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ABSTRACT 


A  sea,  swell  and  surf  program  is  improved,  tested  and  evaluated  on  a 
micro-computer  (HP-9845B).  Sea  swell  is  calculated  by  a  two  dimensional 
spectral  model.  The  energy  balance  equation  is  tested  for  different  cases  of 
wind     velocities     and     water     depths.  Satisfactory     agreement     is     observed 

between  the  offshore  model  and  expected  wave  heights  for  a  15  knot  wind, 
but      the      model     overbuilds      wave      energy     for      a      30      knot      wind.  Wave 

transformation  is  described  by  a  one  dimensional  random  wave  model  in  which 
the  wave  heights  are  described  using  the  Rayleigh  distribution.  The  obtained 
solution  of  the  random  wave  field  is  used  to  predict  the  longshore  currents. 
An  empirical  formula  for  determining  the  breaker  parameters  is  developed, 
based  on  beach  slope  and  incident  wave  steepness.  The  improved  model  is 
tested  using  an  undulated  bathymetry  to  validate  the  model  physics.  The 
model  outputs  of  wave  height  and  current  are  compared  with  data  acquired 
from  a  wave  tank  and  natural  beaches.  The  model  is  found  to  accurately 
forecast  wave  heights,  breaker  location,  breaker  type  and  longshore  currents 
for     several      sets      of     conditions.  Model      limitations      are     discussed      and 
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I.    INTRODUCTION 

A.       OPERATIONAL    REQUIREMENT    FOR    WAVE    FORECASTING 

Wave  and  surf  forecasting  has  been  an  important  part  of  Naval 
Amphibious  operations  since  the  beach  landings  of  World  War  II.  The  need 
for  the  prediction  of  ocean  and  surf  waves  for  military  purposes  led  to  the 
first  real  attempt  to  quantitatively  model  ocean  and  beach  waves.  Ocean 
engineers  and  physical  oceanographers,  such  as  M.P.  O'Brien,  H.U.  Sverdrup 
and  W.H.  Munk,  made  great  strides  in  sea,  swell  and  surf  modeling  during 
the  1940  -  1955  time  period  (Bascom,  1980).  The  wave  and  surf  predictions 
made  by  the  U.S.  Navy  are  still  based  on  these  works.  During  the 
intervening  time  period,  significant  advances  in  knowledge  of  ocean  wave 
dynamics  have   occurred,    and  it   is  time   to   update     prediction   techniques. 

Although  warfare  technology  has  advanced  along  all  fronts,  the  standard 
amphibious  beach  assault  is  still  a  viable  and  preferable  method  of  projecting 
forces  ashore  for  many  tactical  scenarios.  HQwever,  amphibious  assault  craft 
are  vulnerable  to  many  types  of  wave  dynamics.  In  general,  landings  are 
affected  by  the  initial  ocean  sea  and  swell  where  the  boats  are  launched  as 
well  as  by  the  breaking  waves  in  the  surf  zone,  along  with  the  accompanying 
longshore  currents .  Since  in  a  major  landing  tens  of  thousands  of  lives  are 
at  risk,  an  accurate  assessment  of  the  sea,  swell  and  surf  is  vital  to  a 
successful  amphibious  landing  operation    (Joint   Surf  Manual,    1976). 

In  addition  to  the  assault  landing  itself,  an  amphibious  operation  normally 
concludes  with  the  construction  of  a  temporary  harbor  facility  to 
accommodate    the    resupply    of  the    troops  ashore.      The  survivability   and  design 


characteristics  of  these  structures  depend  in  large  part  on  the  environment 
around  them.  The  ability  to  model  wave  height  and  direction,  and  longshore 
currents  is  necessary  to  operational  planners  and  coastal  engineers  in 
designing   and  emplacing  these  facilities. 

In  the  past  twenty  years,  separate  models  have  been  developed  to  deal 
with  ocean  waves,  breaking  surf  and  longshore  currents.  In  1981,  the  U.S. 
Navy  contracted  to  have  a  model  developed  which  would  combine  three 
complementary  models.  The  goal  of  this  comprehensive  Sea,  Swell  and  Surf 
Program  (SSSP)  model  was  to  predict  the  wind  generated  wave  height, 
calculate  the  surf  zone  characteristics  based  on  the  wind  generated  energy 
propagation,  and  use  a  two  dimensional  grid  to  compute  the  longshore  current 
field.  The  model,  developed  by  Wang  and  Chen  (1983),  was  written  in 
FORTRAN  and  designed  to  run  on  a  main  frame  computer.  To  facilitate  the 
model's  use  by  Naval  Oceanographers  in  the  fleet,  the  model  was  converted 
into  BASIC  and  implemented  on  the  HP-9845B/275  micro-computer 
(Devendorf,  1985).  The  primary  goal  of  this  thesis  is  to  modify  the  SSSP 
model  to  improve  its  computational  efficiency. 
B.       WAVE   THEORY    BACKGROUND 

The  water  surface  elevation  in  the  ocean  can  be  viewed  as  a  large 
number  of  individual  wavelets  at  various  frequencies.  The  waves  that  con- 
tribute the  most  energy  to  the  ocean  wave  system  are  wind  generated  waves 
with  periods  of  1  to  30  seconds.  The  size  and  period  of  these  waves  are  a 
function  of  the  velocity  of  the  wind,  the  duration  of  the  storm  and  the 
distance   over   which   the   wind   blows   (fetch). 


The  simplest  theory  is  linear,  small  amplitude  wave  theory.  Sea  water 
is  assumed  to  be  incompressible  and  homogeneous,  and  surface  tension  forces 
are  assumed  to  be  negligible.  The  solutions  of  this  theory  describe  the 
behavior  of  the  individual  components  of  ocean  wind  waves  (Phillips,  1977). 
As  waves  enter  into  shallow  water,  wave  speed  varies  with  the  local  water 
depth  and   refraction  occurs. 

Initially ,  refraction  of  monochromatic  waves  was  modeled  using  wave  ray 
theory    of   a    single    wave    (Arthur,    1949).  It    was   not    until   the    mid    1960's 

that  the  development  of  fast  computers  and  the  implementation  of  better 
physics  enabled  numerical  analysts  to  study  the  generation  of  wind  waves 
over  an  entire  frequency  spectrum.  Early  researchers  in  the  wave  spectral 
analysis  field  introduced  direction  variables  to  study  straight  and  parallel 
contour  wave  spectrum  transformation  (Karlsson,  1969)  and  included  the 
effects  of  bottom  dissipation  and  local  wind  generation. 

As  computers  advanced,  more  complex  refraction  models  were  developed. 
Noda  et  al.  (1974),  using  a  relaxation  finite  differencing  scheme,  solved  for 
the  stationary  wave  spectral  transformation.  The  model  that  is  under  study 
by  this  thesis  was  originally  developed  by  Wang  and  Yang  (1981)  and 
included  the  effects  of  bottom  friction.  Chen  and  Wang  (1983)  improved  the 
model  by  including  the  effects  of  non-stationary  waves. 
C.       CURRENT    MODELING   TECHNIQUES 

Two  numerical  modeling  techniques  are  currently  used  in  wave  and  surf 
forecasting  applications.  The  finite  element  method,  while  highly  complex, 
has  the  advantage  of  being  able  to  incorporate  irregular  localized  bathymetry 
and  has  stable  non-linear  convergence  properties.  Wu  and  Liu  (1985) 
introduced   the   method  in   their  study  of  wave  induced   nearshore   currents. 
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The  finite  difference  method  is  a  widely  used  numerical  method.  It  is 
based  on  a  linear  differencing  approach.  Finite  differencing  normally  uses  a 
rectangular  computational  grid.  This  method  was  used  by  Noda  et  al.  (1974) 
and  Shaiu  and  Wang  (1977)  in  their  studies  of  nearshore  circulation  and  wave 
energy  transformation,  respectively.  The  method  is  easier  to  implement  than 
the  finite  element  method  due  to  its  use  of  a  regular  grid  describing  the 
bathymetry,  and  is  employed  in  the  development  of  the  SSSP  (Wang  and 
Chen,  1983). 
D.      OBJECTIVES 

The  objective  of  this  thesis  is  to  improve  a  sea,  swell  and  surf  program 
which  is  implemented  on  the  HP  9845-B/275  micro-computer.  The  program  is 
divided  into  subroutines,  which  include  the  open  water  wave  generation 
module,  wave  transformation,  the  surf  zone  breaker  calculations  and  the 
longshore  current  routines.  To  test  the  model,  the  output  from  the  improved 
model  is  compared  with  data  bases  acquired  from  two  natural  beaches  as  well 
as  several  wave  tank  experiments.  A  users  manual  has  been  prepared  and  an 
attempt  has  been  made  to  interface  graphic  output  with  the  current  model. 
Limitations  of  the  model  are  discussed  and  future  improvements  suggested. 
This  model  is  proposed  for  use  by  operational  planners  and  Naval  Ocean- 
ographers   at   sea . 
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II.      OVERVIEW   OF  THE   SSSP  MODEL 

A.       OVERALL    PROGRAM    DESIGN 

The  Sea,  Swell  and  Surf  Program  (SSSP)  was  developed  by  Wang  and 
Chen  (1983).  The  SSSP  is  essentially  three  complementary  models  which  are 
merged  to  calculate  offshore  wave  height  fields,  surf  zone  information  and 
longshore  current  velocities.  The  three  modules,  ocean,  surf  and  longshore 
currents,  can  be  run  independently  or  the  offshore  wave  energy  can  be  used 
as  input  to  the  surf  and  current  models.  The  SSSP  is  a  numerical  model 
which   uses  finite  differencing   to  solve   the  governing   equations. 

The  original  version,  as  written  by  Wang  and  Chen,  was  written  in 
FORTRAN  and  designed  to  run  on  a  mainframe  computer.  The  version 
discussed  in  this  thesis  was  converted  to  BASIC  and  implemented  on  an  16 
bit  HP  9845-B/275  micro-computer  (Devendorf,  1985).  The  micro-computer 
version  consists  of  a  control  program  with  subroutines  to  calculate  wave 
direction,  height,  and  number,  and  output  surf,  swell  and  current 
information. 

The  SSSP  considers  an  offshore  wind  generated  energy  spectrum  which  is 
used  to  estimate  the  wave  growth  based  on  initial  inputs  of  wind  speed  and 
wave  height  or  spectral  energy  bands.  As  the  waves  traverse  down  the  user 
defined  grid,  changes  in  the  deep  water  energy  spectrum  due  to  bottom 
friction  and  irregular  bottom  changes,  are  taken  into  account.  If  the  surf 
zone  model  is  to  be  run ,  a  significant  wave  height ,  calculated  from  the  area 
under  the  wave  spectrum,  and  the  peak  frequency  are  used  to  provide  a 
monochromatic    wave    height    input   into    the   surf  model.      It   is   also   possible    to 
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input  monochromatic  wave  data  directly  into  the  surf  module  to  compute  the 
surf  zone  characteristics  without  running  the  computationally  intensive 
offshore   module    (Wang   and    Chen,    1983). 

B.       SSSP    FLOWCHART 

The  SSSP  is  controlled  by  one  main  program  which,  in  turn,  can  call  any 
of  twelve  subroutines.  A  general  flowchart  of  the  entire  SSSP  shows  that 
the  SSSP  is  divided  into  two  main  modules  (Fig.  2.1).  The  user  is  prompted 
for  all  required  information  in  a  menu-driven,  interactive  fashion.  Initially, 
the  user  chooses  whether  to  run  the  offshore  or  the  surf  model.  General  and 
specific  input  is  requested  for  the  module  being  run.  The  open  water  or  surf 
condition  modules  are  calculated  and  the  proper  output  is  generated  by  either 
the    OUTSW   or   OUTSRF   routines. 

A  more  detailed  flowchart  of  the  open  water  condition  module  is 
presented  as  Figure  2.2.  OPEN  computes  the  wavenumber  as  a  function  of 
depth  and  then  loops  through  a  calculation  series  for  each  energy  band.  An 
adjustment  is  made  for  the  angle  of  the  swell/ waves  with  respect  to  the  grid 
axis.  The  DIRECT  and  HEIGH2  subroutines  are  called  and  the  energy  from 
each   spectral  band   is  summed   and  converted   to    wave   height. 

The  original  SURFCON  subroutine  (Figure  2.3)  used  the  most  energetic 
energy  band  computed  by  the  OPEN  module  as  its  frequency  input,  or  used 
data  entered  directly  by  the  user.  The  wavenumber  was  computed,  the 
DIRECT  and  HEIGH2  subroutines  were  called,  the  average  wave  height  was 
calculated   and   a  series   of  breaker   height   and   locations    were   predicted. 

If  longshore  currents  were  desired,  the  NEAR  CI  R  subroutine  was  invoked. 
The    model,    as    implemented    on    the    HP-9845 ,     used    an    iterative    process    on    a 
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Figure  2.2.   Flowcharts  of  the  Open  water,  wave  Direction  and  wave 
He i ant  subroutines  (OPEN,  DIRECT  ano  HEIGH2). 
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Figure  2.3.   Flowchart  of  the  original   Surf  Condition  subroutine 
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two  dimensional  grid   to  solve  for   the  mean   horizontal  velocities.      It    was 
found    that    this   scheme    was    computationally    unstable    as    well    as    requiring    as 
long  as   twenty-three   hours  to  converge   to  a  solution    (Devendorf,    1985). 

To  take  advantage  of  the  most  recent  advances  in  breaker  height 
modeling,  and  increase  the  computational  speed,  the  SURFCON  and  its 
ancillary  subroutines  were  revised  to  use  a  one  dimensional  random  wave 
model.  A'  detailed  flowchart  of  the  new  SURFCON  module  (Fig.  2.4),  which 
replaces  the  original,  shows  the  simplifications  implemented  by  the  new 
model.  DIRECT 2  is  a  simple  Snell's  law  calculation,  applicable  over  straight 
and  parallel  contours.  HEIGH3  uses  a  random  wave  model  which  incorporates 
the  probability  density  function  for  breaking  waves.  The  output  of  the  new 
SURFCON  module  includes  breaker  height,  number  of  breakers,  wave  direction 
and  surf  zone   width. 

LONSHOR  is  a  new  longshore  current  subroutine  (Figure  2.4),  which  was 
implemented  in  place  of  the  unstable  NEAR  CI R  to  simplify  and  speed  the 
calculation  of  the  longshore  currents.  This  model  uses  a  simple  radiation 
stress  balance  to  calculate  the  current  parallel  to  the  beach.  The  next 
section  will  discuss  the  physics  of  the  SSSP  modules  in  detail. 
C.      PROGRAM    MODULES 

1.      Offshore   Module 

When  the  SSSP  is  run,  a  menu  of  initial  choices  is  presented.  If  the 
user  chooses  the  "OPEN  WATER"  option,  the  OPEN  subroutine  is  invoked. 
This  offshore  module  computes  the  wave  height  and  direction  fields  as  a 
function  of  water  depth,  wind  speed  and  direction,  fetch,  and  the  initial 
energy  field.  The  initial  energy  can  be  described  by  a  single  significant 
wave     height,     direction     and     period,      or     as     an     implicitly     entered     energy 
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I  NPUT2 


INPUT3 


rNEWSURFCON 


COMPUTE  WAVE   NUMBER, 
CELERITY,   GROUP  VELOCITY 

CALL  DIRECT2 

CALL  HEIGH3 

CALC  THEORET  I  CA  L   MAX 
BREAKER  HEIGHT 

CALC   FIRST  BREAKER   LINE 
HEIGHT,   LOCATION   AND 
DEPTH 

DETERMINE   BREAKER  TYPES, 
NUMBER  OF  BREAKER   LINES 

CALC   EFFECTIVE   SURF 


I 


OUTSUR  F 


D  I  RECT2 

COMPUTE   WAVE   NUMBER   AND 

DEEP   WATER  CELERITY   C 

o 

DETERMINE   SNELL'S   CONST 
CALC   WAVE   ANGLES 
RETURN   SNELL'S  CONST 
RETURN 


/  HE  I  GH  3  \ 


COMPUTE   WAVE   NUMBER,   C 
AND   Cg 

CALC   INITIAL  ECg ,  £fc 

CALC   I  NTER  I  OR   ECg  ,  £fc 
VALUES 

COMPUTE  RMS   WAVE   HEIGHTS 

CALL   LONSHOR 

RETURN 


LONSHOR 


DEFINE  BOTTOM  SHEAR 
STRESS   COEFFICIENT 

STORE   VELOCUTY  CONSTANT 

CALC  MEAN 
CURRENT  VELOCITY  AT   EACH 
GRID  POINT 

RETURN 


Figure  2.4.   Flowcharts  of  NEWSURFCON,  D1RECT2,  HE1GH3  and  LONSHOR 
subroutines  used  in  the  new  model. 
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spectrum.  Water  depth  for  the  matrix  can  be  read  in  from  a  data  file  or 
entered  directly  from  the  keyboard.  A  tidal  height  correction  constant  can 
be  added  to  the  depth  field  to  take  into  account  the  state  of  the  tide  . 
The  bottom  type  can  be  specified,  which  determines  the  coefficient  of  friction 
for  the  wave  height  calculations. 
a.      Wave   Refraction 

Once  all  the  required  input  is  entered,  the  program  begins  its 
calculation  cycle.  The  first  field  that  is  calculated  is  the  wave  direction 
matrix,    using   a   steady   state   conservation   of  wavenumber    (Phillips,    1977): 

!rr(kcos6)    -    -^(ksine)       =      0  (1) 

o  A  a  i 


where : 


k      =       wavenumber         =      -=—       (m      ) 

±j 

L   =   wavelength    (m) 

6   =   wave  direction  (clockwise  angle  from 
the  X  axis  to  the  wave  ray  (see 
Figure  7) ) 

X   =.  offshore  grid  spacing  (user  specified 
(see  Figure  5) ) 

Y   =   longshore  grid  spacing  (user  specified 
(see  Ficure  5) ) 


The   wavenumber  is   calculated,    for   the   frequency   band   being   considered,    using 
the   linear   dispersion   relation: 

c       =        (g    k    tanh    (k    h))1/2  (2) 
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wnere : 

c   =   angular  frequency   =  2r,i       (sec   ) 

-2 
g   =   acceleration  of  gravity   (m  sec   ) 

h   =   water  depth  (rn) 

f   =   frequency   =   =•  (sec   ) 

T   =   wave  period  (sec) 

Since  (2)  is  transcendental  in  k,  a  sixth  order  polynomial  fit  of  (2)  is  used 
(Hunt,  1979)  to  provide  a  faster  method  of  computing  k  than  the  iterative 
Newton's  method  used  in  the  original  Wang  and  Chen  (1983)  version  of  the 
SSSP.      The    wavenumber   k   is  given   as: 

k      =      ^T7j\-S    +     (1    +    0.66667S    +    0.35550S2 

(gh)X/Z 

+    0.16084S3    +    0.06320S4'-    0.02174SD 

+    0.00654S6    +    0.00171S7    +    0.00029S8 

9-1    V2 

+    0.00011S    )        ]  (3) 


where : 


S      =      o2[|] 

9 


Once  the  wavenumber  is  determined,  subroutine  DIRECT  is 
called  to  compute  the  wave  refraction  of  the  wave  components  over  the 
specified  grid.  DIRECT  starts  with  initial  values  of  the  wave  angles  for  each 
grid    point,    and    uses    a    finite    differencing    scheme    developed    by    Noda    (1972) 
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to  solve  the  differential  equation  for  fl  .  .,  The  differential  equation  is 
center  differenced  in  the  offshore  (X)  direction  and  forward  differenced  in 
the  longshore  (Y)  direction.  The  finite  difference  scheme  has  greater  weight 
on  the  forward  differencing  term  to  increase  computational  stability  (Abbott, 
1979)    and   is   given    as: 

6.     .      =      sm^ir^-Iilksine.    ,  )    -  (1-2t)  (k  sin  6 .  , ,     .) 

l/J  K.  .       .  2.~r  ±  ,  j      1  1  +  1/J 

1  i  J 


+    T(ksin6i+lfj+1)    -   ^(kcosei/j+1) 


(•k  cos  6-  •  ,  )  ]  }  (4) 

1  i   J   -L 


where : 

t   =   weighting  factor   =   0.25 

A  first  guess  for  the  wave  angle  is  calculated  using  Snell's  law 
solution  for  straight  and  parallel  bottom  contours.  The  numeric  scheme  is 
iterated  until  S  .  .  is  within  0.005  of  the  previous  calculation.  The  10  X  13 
grid  point  system  used  by  the  SSSP  is  described  by  Figure  2.5.  The  initial 
differencing  scheme  used  by  the  model  is  a  simple  forward  stepping  scheme 
(Fig.  2.6).  The  convention  for  the  angle  -  is  illustrated  in  Figure  2.7. 
Theta  is  measured  with  respect  to  a  ray  drawn  perpendicular  to  the  beach. 
Therefore,  a  normally  incident  wave  would  have  a  Theta  value  of  180  . 
After  the  direction  matrix  has  been  specified,  the  subroutine  HEIGH2  is 
called  to  calculate  the  height  transformation  as  the  wave  moves  down  the 
grid    (towards  the   shoreline). 
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Figure  2.6.   Numerical  differencing  scheme  used  by  the  OPEN  water 
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Figure  2.7.   Theta  convention  used  by  the  model. 


23 


2.      Wave   Height  Transformation 

This  subroutine  takes  into  account  the  wave  direction  (previously 
calculated),  the  wind  generation,  and  the  bottom  dissipation.  HEIGH2 
employs  the  Phillips- Miles  growth  mechanism  for  the  wind  generation 
calculations.  An  initial  energy  field  is  developed  for  each  grid  point  by 
assuming   that   energy   is   conserved   over  straight   and   parallel   contours: 

£rr(E  Cc   cos    6  )       =       0  (5) 

O  A 


where 


E    =   wave  energy  (joules) 
Cg   =   group  velocity  (m  sec 


The   initial  energy   field   is  derived   by   integrating    (5)    to   yield 


Cg       cos    6 

E      -      E      -12   £  (6 

o     Cg      cos    6  lD 


where   the   subscript    (o)    is  the   initial   condition   at   the   offshore   boundary. 

To    simplify    and    increase    the    model's   efficiency,    Chen    and    Wang    (1983)    make 

four   assumptions: 

a.  Winds   are   steady 

b.  Wave    energy    within    a    spectral    band    is    restricted    to    that 
band      (Allows  superposition   of  wavelets). 


24 


c.  Each     frequency     component     can     be     described     by     a     single 
mean   direction,  G  . 

d.  No    wave-current    interactions. 

These    four    assumptions    allow    the    steady    state    energy    flux    for    each    spectral 
energy   band   to   be    written   as    (Longuet-Higgins   and   Stewart,    1960;    19C2): 


ug    4tt(S  (f  )  Cg  cos  -)     +    yq   £77  (S  (  f )  Cg  sin  6) 

0  A  3  1 


d(f)     +    c,(f)  (7) 


wnere 


S(f)    =   spectral  energy  flux 

_  2 
P       =   density  of  fluid  (kgm   ) 

:,(f)   =   energy  dissipation  due  to 
a  .     - 

DOttorfi  motion 


E  (f)   =   energy  generation  due  to 
wind  stress 

=   mean  wave  direction 


The     wind     generation     term      z        is     the     sum     of    the     linear     and     exponential 
contributions: 


S      =      Sd    +    Sm  <8) 


zo 


wnere : 


S    =   Phillies  Generation   =   oga. 
Sm   =   Miles  generation   =   pg$[S(f)] 


The  energy  generation  term  is  based  on  the  Phillips-Miles  growth  mechanisms. 
Starting  with  a  flat  ocean,  wave  growth  is  initially  linear  until  larger  waves 
are  present,  alter  which  the  growth  is  exponential.  The  linear  growth  term 
(Phillips,  195  7)  is  based  on  random  atmospheric  pressure  variations  acting  on 
a  perturbed  sea  surface.  The  SSSP  first  calculates  the  pressure  fluctuation 
term   as   a  function   of   wavenumber   and   frequency: 

P(k.  „).    6-13.'""^6;^ '1 Tv  vi 


2J  '  2  2 

+  (k  sine)    v  .  +  (  k  cos  :  -A ) 


(9 


wnere 


W    =   wind  soeed  in  knots 


y,   =   0. 


W 

,1.28 


T 


G  .  95 
v 2   =   0.52  A 

=   angle  between  wind  and  wave 

After   P(k,  0  )    is   calculated,    the   linear  growth   term,  Q  ,    is  determined   by 


4tt    kc 
a       =      2~T~    P(k'    c  )  dO) 

p  g 

w^3 
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where : 

_  -> 
p    =   density  of  water  (kg  m      ) 

The  exponential  growth  term  becomes  important  after  waves  have  formed 
(LaBlond  and  Mysak,  1978).  When  flow  separation  occurs  in  the  lee  of  a 
wave  crest ,  momentum  is  transferred  to  the  wave  because  of  the  pressure 
differential  resulting  in  exponential  wave  growth  (Miles,  1957;  1959a, b; 
1962) : 


(^-)  [(|)COS  6  -b]  (11; 


where 


a   =   constant   =   5.0 


S   =   ratio  of  air  to  water  density 
C   =   wave  celerity  (m  sec   ) 
b   =   constant   =   0.90 


The  Phillips-Miles  wind  wave  generation  has  no  mechanism  for 
shutting  off  the  wave  growth.  In  nature,  wind  waves  will  grow  as  a 
function  of  wind  speed,  duration  and  fetch.  But  after  growing  to  a  certain 
size  and  steepness,  the  waves  will  become  unstable  and  break.  Therefore,  it 
is  necessary  to  have  a  limiting  function  so  that  the  waves  will  only  grow  to 
a  "fully  arisen"  sea.  Several  functions  have  been  proposed  to  introduce  an 
"equilibrium  value",  above  which  the  wave  growth  stops.  The  SSSP  uses  a 
Pierson  and  Moskowitz  (1964)  fully  arisen  sea  in  deep  water.  The  high 
frequency      tail     is     described      by      a      Phillips      (1957)      equilibrium     condition, 


modified   by   Kitaigorodskii   et   al.     (1975)    and   Thornton    (1977)    for   intermediate 
or  shallow   water: 


S(f , 6)        =       B  g2  f       H(o)  r(W    )  (12 


where : 


8      4 
H(o)    =   j  7i  cos  6   accounts  for  angular 

spreading 
_9  2  W   I (W  )    -1 

nwn)     =    :  2(wr)  [i  +  L— n — ] 

sinh[2W  I (W  ) ] 

n    n 

Wn     =   0(h)  V2 
n         a 


I(W    )    meets   the   condition   that      I   tanh(W     "  )    =   1 
n  n  ~ 


2 
7  (W    )   =   1    in   deep    water    and   for    W      <    1,    ?(W    )  —  1/2    W        and    S(f,  c    )   = 
n  n  n  n 

-3  ~9 

1/2    q   g   h   df        H(  6  )    where  q   is   on   the   order   of  10    ~    (SSSP   uses   a   value   of 


0.073). 


Bottom  dissipation,  z  ,  is  calculated  as  the  work  done  on  the 
bottom  due  to  the  bed  shear  stress.  Assuming  the  usual  quadratic  bed  shear 
stress    law,     the    wave    energy    dissipation    for    a    particular    spectral    frequency 
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band   is  described   by   Hasselmann   and    Collins    (1968) 


£d(f)   =   pwCf  l°bl  Ub(f) 


where : 

Cr      =   bed  shear  stress  coefficient   =   0.01 

I U,     =   total  flow  field  bottom  velocity 

U,  (f)   =   wave  induced  water  particle  bottom 
velocity  for  a  specific  spectral 
component 

Using  linear  wave  theory  transfer  functions  to  relate  the  surface  elevation 
energy  spectrum  to  the  bottom  velocity  field,  the  frequency  dependent 
dissipation  function      tAt)    can   be   expressed   as: 


Ed(f)       =       P  Cf  T(f)    S(f)    |Ub|  (14) 

111  2 

where   the   transfer  function     T(f)    =    [    g     k   /  c      cosh      (kh)].    The   total    wave 

induced     bottom     velocity      |  L\    |     is     calculated     by     integrating     the     spectral 

transfer  function   across   the   entire  frequency   spectrum: 


1/2 
Ub!       =       I  /     T(f)    S(f)  df]  C5) 

0 


Bottom  dissipation  in  the  SSSP  is  only  calculated  when  deaLng  with 
intermediate  or  shallow  water  waves.  If  deep  water  wave  heights  are  being 
computed,    the   value     z      is  set   equal   to   zero    (Wang   and    Chen,    1983). 
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3.  Longshore    Currents 

One  of  the  most  serious  problems  with  the  original  model  was  the 
instability  of  the  nearshore  circulation  calculations.  As  much  as  23  hours  of 
CPU  time  were  required  for  the  subroutine  to  converge  to  an  answer.  It 
was  decided  that  a  one  dimensional  longshore  current  model  which  assumes 
straight  and  parallel  bottom  contours  would  be  substituted  for  the  unstable 
two-dimensional  model.  This  model  uses  a  simple  alongshore  radiation  stress 
balance  to  calculate  the  mean  longshore  current.  Although  no  circulation 
information  is  derived,  the  new  longshore  current  subroutine,  LONSHOR, 
provides  the  longshore  current  velocity  efficiently  and  accurately.  The 
physics   and   details  of  the   LONSHOR   subroutine   are  outlined  in    Chapter   II. 

4.  Surf  Zone 

The  original  SURF  CON  module  made  calls  to  the  same  wave  direction 
(DIRECT)  and  wave  height  (HEIGH2)  subroutines  as  the  offshore  module. 
One  disadvantage  of  this  scheme  is  that  these  subroutines  are  very 
computationally  intensive,  especially  when  considering  the  narrow  width 
between  grid  points.      A   faster   wave   height   model   was   needed. 

The  state  of  the  art  in  breaker  height  models  uses  a  probability 
density  function  to  describe  breaking  wave  heights  (Thornton  and  Guza, 
1983).  Some  limitations  are  imposed  by  using  this  new  model,  however. 
The  model  assumes  that  the  beach  has  a  local  straight  and  parallel  bottom 
contour.  This  bathy  metric  restriction  allows  longshore  features  such  as 
sandbars  to  be  modeled  accurately,  but  offshore  feature  such  as  channels  or 
canyons  are  smoothed  out  by  a  depth  averaging  subroutine.  It  is  felt  that 
this    assumption    is    reasonable,     in    light    of    the    fact     that    many    beaches    of 
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operational    interest    are    of   the    "straight    and    parallel"    type    within    the    first 
two   hundred   meters  from   shore. 

The  above  assumption  means  that  a  simple  Snell's  law  calculation  of 
the  wave  direction  is  possible.  A  short  subroutine,  called  DIRECT  2,  is 
implemented  to  calculate  the  wave  direction  matrix.  The  use  of  the  new 
height  and  direction  modules,  HEIGH3  and  DIRECT2  have  increased  the 
efficiency  of  the  SURFCON  calculations  immensely.  The  average  run  time  for 
the  surf  zone  model  predicted  wave  heights  has  decreased  from  over  four 
minutes  to   less   than   30   seconds. 
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HI.      MODEL   IMPROVEMENTS 

A.       OFFSHORE    WAVE    HEIGHT    GROWTH 

The  current  version  of  the  SSSP,  converted  to  run  in  HP  BASIC  on  the 
HP-9845B/275  micro-computer,  gives  reasonable  values  for  wave  height 
growth  for  wind  speeds  of  15  knots  or  less  (Devendorf,  1985).  However,  for 
higher  wind  speeds  of  30  knots,  the  SSSP  exhibits  a  growth  rate  that  is 
significantly  faster  than  the  growth  of  the  JONSWAP  spectrum  (Bishop, 
1983),  used  for  comparison  in  this  analysis.  The  JONSWAP  spectrum  was 
chosen  as  a  basis  for  comparison  because  the  JONSWAP  curves  were 
experimentally  derived  from  the  limited  fetch  North  Sea  area.  It  is  foreseen 
that   the   SSSP   model  may   be  applied  to  similar  operational  areas. 

For  30  knot  winds,  the  model's  wave  growth  not  only  occurs  too 
quickly,  reaching  a  maximum  value  within  50  nautical  miles,  but  the 
magnitude  of  the  growth  is  too  high  (maximum  wave  height  of  6.0  meters 
compared  to  a  JONSWAP  maximum  of  4.2  meters.  Before  sensitivity  tests 
were  run  on  the  wind  generation  portion  of  the  model,  the  code  was  checked 
carefully    to    ensure    that    the    equations    were    being    implemented    correctly. 

To  determine  the  sensitivity  of  the  model  wave  growth,  several  runs 
were  made  with  changes  to  the  key  growth  parameters.  These  parameters 
include  A,  the  Phillips  linear  growth  term,  Equation  9  (Phillips,  1957);  B, 
the  Miles  exponential  growth  term,  Equation  10  (Miles,  1957;  1959a, b; 
1962),    and  R,    a   weighting  factor  set   equal  to  1.0  in  the  original  model. 
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The    first     run     was     used     as     a     standard     to     compare     the     modifications 
against.      The  inputs  for  the   OPEN   module   were: 

Grid   Spacing  X      =        Y      =      10    NM 

Wave   Direction        Theta  =  180      (normal  incidence) 

Wind    (W)      =     30   and   15    kts  @   270° 

Shoreline  =  270 

Bandwidth    (Df)      =  0.01    Hz 

Frequency  bands  =0.05  -  0.14  Hz 
The  output  for  these  "standard"  model  runs  are  shown  as  Figures  3.1  and 
3.2.  The  overbuilding  of  the  open  water  waves  for  W  =  30  kts  is  evident 
when  compared  with  the  JONS  WAP  spectrum  wave  growth  (Fig.  3.1).  For  a 
more  moderate  case  where  W  =  15  kts,  the  model  wave  growth  shows 
satisfactory  agreement  with  the  JONSWAP  growth  (Figure  3.2).  To  see  the 
impact  of  changing  the  weighting  factor,  R,  two  runs  were  made  setting  R 
equal  to  0.5  and  0.1  respectively.  The  predicted  wave  height  with  R  equal 
to  1,  0.5  and  0.1.  is  shown  in  Figure  3.3.  Setting  R  equal  to  0.1  (chain 
dashed  line)  gives  a  good  fit  for  high  wind  speeds  but  underbuilds  the  wave 
field    for     lighter     winds     (Fig.     3.4).  Clearly,     a     simple     change     of     the 

weighting  factor   R    will  not   resolve   the  problem. 

The  Miles  exponential  growth  term,  B,  is  in  the  denominator  of  the 
energy  calculation  so  its  value  was  increased  to  determine  the  growth 
retarding  effect  it  would  have.  Two  runs  were  made  with  B  doubled  (chain 
dotted  line)  and  with  B  increased  by  a  factor  of  10  (not  shown)  ,  depicted 
by  Figure  3.5  The  result  of  increasing  the  exponential  growth  term  is  to 
flatten  out  the  upper  portion  of  the  growth  curve  which  is  not  the  correction 
that   needs  to  be  applied. 
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Because  the  model's  rapid  growth  appears  to  occur  in  the  linear  region 
of  the  growth  curve ,  the  Phillips  linear  term ,  A ,  was  analyzed  closely  to 
see  if  a  simple  constant  multiple  would  bring  the  wave  growth  in  line  with 
the  JONS  WAP  curves.  A  constant  multiple  of  0.125  is  applied  to  the  A  term 
(Fig.  3.6).  This  method  gives  a  better  fit  to  the  JONSWAP  curve  for  high 
wind  speeds,  although  the  SSSP  still  over  builds  the  wave  heights.  To  test 
this  multiple  at  lower  wind  speeds,  a  model  run  was  made  for  a  wind  speed 
of  15  knots,  for  a  bandwidth  of  0.05  Hz  and  frequencies  of  0.05,  0.10, 
0.15,  and  0.20.  The  results  of  the  model  run  are  compared  with  the 
equivalent  JONSWAP  curve  and  the  curve  from  the  existing  model  (Fig.  3.7). 
The  fit,  while  not  exact,  is  reasonable,  suggesting  that  the  A  term  is 
responsible  for  the  overbuilding  of  the  model  wave  heights. 

The  actual  cause  of  the  overbuilding  of  wind  waves  by  the  SSSP  is 
unknown.  Several  ad  hoc  "fixes"  to  the  model  are  suggested  with  the 
recommendation  that  further  research  on  this  problem  be  conducted.  A 
simple  constant  multiple  of  the  linear  growth  term,  on  the  order  of  0.10  to 
0.15  is  one  recommended  empirically  derived  modification  to  the  wind 
generation     calculation.  Another     possible     modification     is     to     have     the 

weighting  function,  R,  decrease  as  the  wind  increases,  to  act  as  a  damper 
on  the  over  generation.  At  present,  the  SSSP  wave  generation  routine 
remains  in  its  original  form.  It  is  felt  that  the  model  is  adequate  for  low 
wind  speed  generation.  Operational  planners  should  be  cautioned  about  the 
high  wind  speed  generation  problems. 
B.      WAVE    REFRACTION   MODULE 

To  test  the  SSSP  model's  stability  during  the  wave  height  prediction  in 
the    coastal    zone,    a   sinusoidally   varying   bottom    was    explicitly   entered   using 

38 


5t} 


>.> 


&30 
>  r5 
.<  ii 


£ 


E=;  ii 

g 

>E 

00 

C 


3c  Oo  8T  9T  tT  2T  OT    8     9     t     2 

(li)  1HOI3H  3AYAV 


r    ^c 


— *    i 


OT 

— 


re 
— 


bo  L    O) 

C      -5      L. 


«     OT 

C 

ii      0) 

L. 


«  Z- 


c    c 


oE 


C    C/2 


"    2 

.  o 


OT     t>C 

—    _ 

re 

>>: 

>    a) 


.So- 


ot — 

C  x 

o  o 

CO 


OT 


0) 

u 

a 


39 


§ 


r-<0 

2  H 
'£< 
to 

>> 

<Joo 

o 

"-3 


S3  02  8T  91  *T  Z\  0T    8     9     t     Z 

CldO  1HDI3H  3AVM 


40 


ss 

So 

s » 

D-  CO 

cow 

C0>< 

cob 
>> 

<Jco 

^2 

O 

•-5 




7 

■■■■(}■- 

■x 

— -o» 

....A,.. 

£ 

r 

I  | 
1  | 
1  ' 
1  ' 
1   ; 

m, i 

i 
i  • 

i  \ 

1   • 
1   ; 

1   ' 

K 

h 

:  \ 
:    \ 

\ 

\ 

\ 

\ 

\ 

'.             \ 
.  :            \ 
V             \ 

V, i 

\ 

\ 

|u 

\ 

\ 

f 

\ 
\ 
\ 
s 
\      \ 

\\ 

•                             1 

\ 

v\ 

i 

\ " 

! 

f^ 


3 

52 


O 


r-i 

u 

< 


CO 

5  5 


r- 


o 


''■  ~z 


r. 

Ifl     0) 
—  *0 


_    9, 

~  o 


.=      ~ 


o  -c 


o 


~   o 
o  E 

ii   a. 

en 
<  «; 

tn 


m  X  Z 

0>  «  S 

«->  be  c 

>>  £  w 

r  a  c 

— <  cc 

09 

e  c  a) 

.  —  _ 


0) 

3 

be 


g       t       e       z       i       o 
(II)  1HDI3H  3AYM. 


41 


the  following   equation: 

h,  v        =        0.025X    +     Eh.    (r~)  (1    "     (rr-))COS(^)  (16) 

(x,y)  b    Xb  Xb  Ly 


where 


£    =   bottom  perturbation  term  (0.2) 

h,   =   breaker  depth  (2.5  m) 
b  r 

X,   =   offshore  distance  to  breaker  line 
(100  m) 

L    =   wave  length  of  bottom  perturbation  (240  m ) 

This  equation  provides  a  smoothly  varying  bottom  profile  in  both  the  X 
(offshore)  and  Y  (longshore)  directions.  To  check  that  the  bottom  profile 
was  being  accepted  correctly  by  the  model,  and  that  the  calculated  results 
were  reasonable,  approximately  25  model  runs  were  completed  with  varying 
input  parameters.  Selected  cases  are  presented  herein  for  discussion.  In 
the  cases  presented,  the  offshore  waves  are  normally  incident  to  the 
shoreline.  Obliquely  incident  model  runs  were  made  with  Theta  ranging  from 
+10  to  -10  degrees  from  normal,  for  comparison  with  the  normally  incident 
case . 

1.       Model  Experiments 

The  input  parameters  in  the  model  are:  significant  wave  height,  H 
(1  m);  wave  period,  T  (12  s) ;  offshore  boundary  depth,  h(l2  m);  artificial 
viscosity,  T  (0.0);  and  bottom  perturbation  term,  £  (0.2).  The  beach 
faces  270  and  the  normally  incident  wave  rays  come  from  270  .  This  gives 
a  Theta  value  of  180  according  to  the  model  convention  shown  in  Figure 
2.7. 
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The  model  results  for  wave  heights  and  directions  are  shown  in 
Figures  3.8-3.15.  The  longshore  depth  profiles  for  locations  inside  (X=3j 
and  outside  (X=8)  the  surf  zone  are  chosen  to  examine  the  wave  refraction. 
Grid  point  X  =  3  implies  a  distance  of  40  meters  offshore,  well  inside  the 
breaker  line,  while  X  =  8  is  140  meters  offshore  and  outside  the  breaker 
line.  The  depth  profile  presented  in  Figure  3.8  shows  that  the  depth  varies 
sinusoidally  in  the  longshore  direction .  The  beach  profile  along  grid  point  Y 
=  3  is  shown  as  Figure  3.9.  As  shown,  the  depths  are  slightly  deeper  than 
the  plane  beach  inside  100  meters  (breaker  position),  and  shallower  outside 
100    meters. 

As  waves  enter  the  surf  zone,  their  increase  in  height  is  a  function 
of  water  depth.  For  the  points  outside  the  breakers,  the  largest  waves  are 
expected  to  occur  where  the  water  depth  is  the  shallowest  due  to  shoaling. 
As  can  be  seen  from  Figure  3.10,  the  peak  waves  (dashed)  are  towards  the 
edges  of  the  domain,  but  they  do  not  coincide  with  the  minimum  in  depths  at 
0  and  240  meters.  This  anomaly  is  a  result  of  the  boundary  conditions  of 
the  model.  Also  note  the  asymmetric  peaks  at  40  and  200  meters.  This 
implies  errors  in  wave  height  prediction  due,  in  part,  to  the  wave  direction 
calculations.      These   errors   will  be   addressed   later  in   this  discussion. 

Wave  direction  is  calculated  in  terms  of  Theta .  An  angle  greater 
than  180  tends  to  turn  the  wave  ray  to  the  left  while  Theta  less  than  18  0 
turns  the  ray  to  the  right  (See  arrows  in  Figure  3.11).  The  model  at  least 
quantitatively  predicts  the  wave  refraction  pattern  correctly  for  all  cases, 
consistently  turning  the  rays  to  "attack"  the  shallower  water  depths  outside 
the  breaker  line.  The  model's  boundary  conditions  force  0  =  8  _  and  Q 
=     9        where     G    is   the    wave    direction    at    the    n        long   shore   grid    point .    This 
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accounts  for  the  departures  of  the  wave  direction  at  the  edges  of  the 
domain.  The  dotted  lines  are  the  interpretation  of  how  the  wave  direction 
patterns  should  behave.  In  the  middle  of  the  domain,  where  both  cases 
exhibit  a  nodal  point,  the  wave  direction  passes  through  180  or  normal  to 
the  beach  as  expected. 

The  final  test  conducted  compared  the  effect  of  Tau  on  the  wave 
direction  calculations.  Tau  is  an  artificial  viscosity  parameter  which  is 
intended  to  be  a  smoothing  term  for  the  calculation  of  wave  angle. 
Changing  Tau  from  0.00  to  0.25  had  a  negligible  effect  on  the  wave  height 
calculations,  but  the  difference  in  direction  increased  with  decreasing 
distance  from  shore.  The  case  where  X  =  8  is  plotted  in  Figure  3.12.  The 
introduction  of  the  Tau  term  (dashed)  decreases  the  magnitude  of  the  wave 
angle,  but  the  maximum  (minimum)  values  are  not  exactly  in  phase  with  the 
Tau  =  0.00  case.  The  maximums  are  shifted  slightly  toward  the  edges  of  the 
domain.      This  is  the  expected  behavior  of  the  smoothing  term. 

The  dissipative  model  exhibits  unexplained  behavior  inside  the 
breaker  line  as  depicted  by  Figures  3.13-3.15.  At  X  =  5  (Fig.  3.13), 
the  Tau  model  reduces  the  magnitude  of  the  wave  direction  field.  At  X  =  4 
(Fig.  3.14),  spurious  asymmetry  appears  in  the  direction  calculations  with 
the  Tau  term  applied.  Furthermore,  at  X  =  3  (Fig.  3.15),  the  effect  of  Tau 
is  to  increase  the  magnitude  of  the  wave  direction.  This  apparent  angular 
deviation  has  no  effect  on  the  model's  wave  height  calculations  because, 
inside  the  breaker  line,  wave  height  becomes  a  function  of  depth  only. 
However,  this  instability  does  effect  the  longshore  current  calculations, 
which  are  sensitive  to  local  wave  direction  inside  the  surf  zone . 
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Another  indication  of  a  numeric  error  in  the  wave  direction 
calculation  was  noted  in  the  exact  center  (Y  =  7)  of  the  domain.  Since  the 
depth  profiles  exhibit  a  maximum  (minimum)  value  at  the  domain  center, 
theory  suggests  that  the  waves  remain  normally  incident  at  these  nodal 
points.  Therefore,  the  Theta  value  should  be  180  degrees  at  the  center  of 
the  domain.  The  model  output,  both  with  and  without  Tau,  exhibits  a 
monotonically  increasing  error  as  seen  in  Figure  3.16. 
2.      Improvement  of  the   Numerical   Scheme 

In  an  effort  to  eliminate  the  numerical  instability  in  the  wave 
direction  module ,  the  finite  differencing  scheme  was  analyzed  to  see  if  an 
improvement  could  be  made.  The  original  scheme,  pictured  in  Figure  2.6, 
uses  a  simple  centered  and  forward  differencing  method  and  is  fairly  robust  in 
the  longshore  ( Y )  direction  but  is  less  so  in  the  offshore  ( X )  direction .  To 
increase  the  accuracy  of  the  numerical  scheme,  the  method  pictured  in 
Figure  3.17  was  substituted  in  the  DIRECT  subroutine  for  the  original 
scheme.  This  particular  method  is  a  well  accepted  technique  which  has 
proven   to  be  stable   and  free   of  damping    (Haltiner  and   Williams,    1980). 

Several  model  runs  were  made  to  test  the  new  technique.  Using  the 
standard  input  parameters  and  setting  Tau  equal  to  0,  the  first  result  was  a 
gain  in  accuracy  of  the  wave  direction  calculations.  The  domain's  central 
directions  were  exactly  180  (to  three  significant  figures),  as  seen  in  Figure 
3.18,  and  the  refraction  angles  were  symmetric  about  the  center  of  the 
domain  as  expected  due  to  the  symmetrically  varying  bathymetry  (Figure 
3.19). 

A  second  result  of  the  implementation  of  the  improved  numerical 
technique    is    the    stabilizing    effect    on    the    wave    height     calculations.        The 
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previous   asymmetric    wave    heights    are   now   symmetric    with    respect    to   height 
and  peak  location    (Fig.    3.20). 
3.      Boundary   Conditions 

The  SSSP  model  imposes  no  flux  boundary  conditions  on  the 
refraction  and  wave  height  calculations  so  that  the  left  and  right  boundaries 
are  specified  as  follows: 

Q(M,1)   =   Q(M,2) 
Q(M,N)   =   Q(M,N-1) 
where:      Q     =  the  quantity   being  calculated 
M     =  offshore  row  number    (10   to   l) 
N     =   the   last   longshore   column    (13) 

The  offshore  row  is  completely  specified  by  the  initial  conditions.  The 
calculations  then  proceed  from  offshore  to  onshore  and  from  left  to  right  on 
the  grid. 

No  flux  boundary  conditions  account  for  the  unusual  behavior  of  the 
wave  height  and  direction  curves  at  the  edges  of  the  domain.  To  satisfy  the 
condition  of  well-posedness  for  the  applied  model,  cyclic  boundary  conditions 
are  substituted  such   that: 

Q(M,1)   =   Q(M,N-1)  (17) 

Q(M,N)  =  Q(M,2) 
A  plot  of  wave  height  versus  longshore  distance  at  X  =  8  using  the  same 
inputs  as  the  previous  case  is  presented  in  Figure  3.21.  Note  the  shift  of 
the  peaks  toward  the  edges  of  the  domain  as  well  as  the  abrupt  increase  in 
the  magnitude  of  the  wave  height.  These  numerical  results  now  match  quite 
well     with     theory     which     requires     that     the     wave     rays     turn     toward     the 

58 


shallowest    water,    implying   that    the    largest    wave   heights   should   occur   at   the 
edges  of  the  domain. 

The     wave     direction     calculations     also     show     improvement.         Close 
agreement    between    the    cyclic    boundary    condition    values    and    the    expected 
values    (dotted)    is  shown  in   Figure   3.22.      Additionally,    the   average   value   on 
the   boundary  segment   agrees   closely   with   the  exact   value    (solid  circles). 
4.      Results 

Based  on  the  results  of  many  model  experiments,  the  new  model 
appears  to  offer  an  increase  in  stability  and  computational  accuracy.  The 
improved  wave  direction  subroutine  consistently  turns  the  wave  rays  in  the 
correct  direction.  The  magnitude  of  the  direction  changes  is  qualitatively 
correct.  The  wave  height  subroutine  works  well  both  inside  and  outside  the 
surf  zone.      Wave  heights  are  highest  at   the  shoal  and  lowest  at   the   trench. 

One  drawback  to  the  original  wave  height  subroutine  is  the 
computational  speed.  When  running  the  program  to  generate  an  energy 
spectrum  for  sea  and  swell,  the  model  takes  approximately  45  minutes  per 
frequency  band.  Even  the  monochromatic  surf  module  is  slow  (4.5  -  5 
minutes)  because  of  the  number  of  iterations  required  to  converge  to  the 
preassigned   criterion. 

Least  acceptable  of  all  is  the  Nearshore  Circulation  module  which 
proved  to  be  unstable  as  well  as  requiring  as  much  as  23  hours  to 
asymptotically  approach  an  answer.  To  provide  fast,  accurate  longshore 
current  information,  a  one  dimensional  model  is  substituted  for  the  two 
dimensional  circulation  model.  While  some  potential  rip  current  information 
will  not  be   provided  by   this  simple   model,    the   one   dimensional    model   has    the 
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advantage    of    being    fast,     accurate    and    able     to    provide     the     most     important 
input,    the   longshore   current   distribution. 

A    state    of    the    art,     one-dimensional     wave    height    model    is    used    as 
input    to    the    longshore    current    model.       The    new    wave    height    and    longshore 
current   subroutines   are   described   below. 
C.       ONE    DIMENSIONAL    WAVE    HEIGHT    CALCULATION 

To  increase  the  SSSP's  computational  efficiency  while  maintaining 
reasonable  accuracy,  a  one  dimensional  random  wave  height  model  was 
substituted  for  the  original  two  dimensional  monochromatic  model.  This 
model  characterizes  the  transformation  of  the  wave  height  probability  density 
function  (pdf)  from  offshore  to  the  shoreline  (Thornton  and  Guza,  1983). 
The  model  assumes  that  the  waves  are  narrow  banded  in  frequency  and 
direction  such  that  the  waves  are  specified  by  a  single  mean  frequency  f  and 
mean  direction  6  ,  that  the  bottom  contours  are  straight  and  parallel,  and 
that  the  wave  conditions  are  stationary.  This  probabilistic  approach  to  the 
wave  breaking  problem  was  originally  suggested  by  Collins  (1970)  and  Battjes 
(1972).  The  form  that  is  used  here  is  based  on  integrating  the  energy  flux 
balance  equation.  The  calculated  wave  heights  are  a  function  of  the  shoaling 
conditions  of  the  integral  path  from  offshore  to  the  shoreline  (Battjes  and 
Janssen,    1978). 

For  straight  and  parallel  contours,  the  energy  flux  is  balanced  by  wave 
breaking   and   the  frictional   dissipation: 

d  (E  Cg     cos    6)  .  (is) 


'6X 


<v  +  <cf> 
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where : 

<£,>   =   dissipation  due  to  wave  breaking 

<£f>   =   frictional  dissipation 

In  the  domain  of  interest,  it  has  been  shown  that  the  frictional  dissipation  is 
less  than  3%  of  the  dissipation  due  to  breaking  waves  (Thornton  and  Guza, 
1983)    and   will  be   neglected. 

The    waves   in    the   new    model    are    described    by    the    Rayleigh    wave    height 
distribution : 


ou      H    '  (19) 

,„,       2H       rms 
p(H)   = e 

rms 


where : 


H      =   root  mean  sauare  wave  height 

rms  -  3 


where    H  is    the    root    mean    square    wave    height .       This    distribution    takes 

rms 

into  account  the  fact  that  there  is  a  distribution  of  wave  heights  throughout 
the  entire  field  for  a  given  frequency.  Previous  random  wave  models 
truncated  the  Rayleigh  distribution  in  various  ways.  Recent  observations 
have  found  that  even  after  the  waves  break,  the  Rayleigh  distribution-  is  still 
valid  (Thornton  and  Guza,  1983).  Therefore,  in  the  improved  model,  the 
waves  are  everywhere  Rayleigh  distributed.  The  probability  of  a  wave 
breaking  Ph(H)  is  described  by  a  simple  weighting  of  the  Rayleigh 
distribution: 

Pb(H)   =   W(H)    p(H) 
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The    weighting    function,     W(H),     which    makes    larger    waves    at    a    given    depth 
more   likely   to   break,    is   given   by    (Thornton    and    Guza ,    1983): 


H  4 

W(H)      =       [-SS*]  (20J 


where 


Y   =   adjustable  parameter  based  on 
field  data 

The  model  assumes  that  breaking  waves  are  similar  to  periodic  bores 
(LeMehaute',  1962)  as  depicted  in  Figure  3.23.  The  bore  dissipation  per 
unit   area  is   calculated: 

1       „     (h2~hl)3.  1       a     (BID3    n 

£bore      "      TPg         05-      Q      "      4pg— 2~Q  (21) 

l    i  n 

where: 


Q   =   volume  discharge  across  bore 

B   =   breaker  coefficient  (described 
from  data) 


The     volume     discharge     parameter     is     described     for     a     linear     periodic     bore 
(Hwang  and   Divoky ,    1970)    by: 

Q   =   Ch/L  (22) 

where    C    is    the    wave    celerity.       Using    linear    wave    theory,    the    wave    energy 
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and   group  velocities   are   given   by: 


00 


E       =      pg    /       H      p(H)   dH       =       i  0  g    h2 


o  rms  '23J 


and, 


C  =      §[1+  2   k  h    ,     cos    -  .24) 


After  integration    and    the   substitution    of  the    probability    density   function, 
Equation   21   leads   to: 


<£,>       =      ^    ^og    -=#--   H  (2a) 

b  16  y      4,d      rms 

Y    h 


where: 

f   =   mean  frequency 

This  equation  for  the  bore  dissipation  can  be  tuned  to  the  data  by  varying 
the  breaker  coefficient  B ,  and  the  adjustable  parameter  )  .  Substituting  the 
bore   dissipation  function    (25),    the   energy  flux    equation    (23)    is   given   by: 


§x<lP*>4sCgcose,     -    ±,-,9J^nlms         Ml 


For  shallow   water,    the   group  velocity    Cg   is  approximately   equal   to   the    wave 

1/2 
celerity    C    which    equals    (gh)         .       Equation    (26)    was    chosen    to    be    used    in 

the     wave    height    module    because,     although    a    more     complicated    expression 
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gave  a  slightly  better  fit  to  the  data,  Equation  (26)  provides  an  analytical 
solution     which     allows     a     detailed     study     of     the     model's     behavior.  The 

analytical  solution  is  obtained  for  normally  incident  waves  impinging  on  a 
plane  sloping  beach  and  is  used  to  compare  the  accuracy  of  the  numerical 
model.  The  analytical  solution  is  solved  by  integrating  (26)  for  the  wave 
height    (Thornton   and    Guza ,    1983): 

y 

H      _    1/5  9/10      h23/4    1       a  vl-l/5 
Hrms   "   a    h     [1  -  h     ^2374  '  ^5/2 )  ]      (27) 

o       o 

where : 

a  lA(Z)1/2     V4  tan  B 

15  *        B3  f 

tan  8   =   slope  of  beach 

y0    =  h2  h1/2 
o         o  o 

The  numerical  model  uses  a  modified  Euler  integration  scheme  to 
integrate  the  energy  flux  equation.  A  weighted  predictor/  corrector  loop  is 
used  to  increase  accuracy  and  stability  as  the  water  depth  becomes  shallow. 
The   basic  calculation  is  a  forward  stepping  scheme   with  a  variable  grid  size: 


EC^x    2      =      EC^x 


1    2    b 


2  b 


)  AX 
1  (28) 


Well  outside  the  surf  zone  the  wave  height  changes  are  very  slight  with 
distance  as  the  bore  dissipation  is  essentially  zero,  and  are  relatively 
insensitive  to  the  grid  size.  However,  once  the  waves  start  to  break,  the 
bore    dissipation    increases    rapidly,    modifying    the    wave    heights.      To   increase 
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the  resolution ,  the  mesh  size  decreases  by  a  factor  of  four  when  the  bore 
dissipation  reaches  a  value  of  0.01.  This  flag  value  was  arrived  at  by  a 
series  of  model  tests  under  different   initial  conditions. 

The    model    was    tested    against    the    analytical    solution    for    a    variety    of 

initial   values.      The  initial  conditions  insured   that  shallow   water  criteria   were 

1/2 
met    so    that    the    shallow    water    approximation    of    Cg    =    (gh)         ,    rather    than 


the    explicit    Equation    (24),    could    be    used.        The    first    test    used    a    shallo 


v\ 


steepness     offshore     wave     (H    /L       =     0.01)     having     a     height     of     1.0     meter 

normally   impinging  on  a  plane  beach   with  a  slope  of  0.04.      The   period  of  the 

wave    was    10    seconds.         The     numerical    and     the     analytical     solutions     are 

compared    in    Figure    3.24.        Agreement    is    very    close    through    most    of    the 

offshore   grid  points.      As  the    water  depth   gets  very  shallow,    the  step  size   is 

again  reduced  to  insure  the  solution  converges  for  the  last  few  grid  points. 

For   a   steep   wave    case    (H    /L      =   0.02),    an   incident    wave   of   2.5    meters 

o     o 

was  used  as  the  offshore  input  (Fig.  3.24).  The  wave  period,  beach  slope 
and  normal  incidence  all  remained  the  same.  The  numeric  and  analytic 
solutions  are  again  in  good  agreement ,  although  the  numeric  solution  slightly 
underpredicts  the  wave  height  as  it  did  in  the  case  where  H  equals  1.0  m. 
It  is  of  interest  to  note  that  the  high  wave  steepness  case  shows  no  growth 
in  wave  height  due  to  shoaling  prior  to  breaking.  The  wave  heights  decay 
monotonically   toward  the  shoreline. 

Having  satisfactorily  compared  the  model  results  with  the  analytic 
solution,  further  improvements  to  the  model's  computational  efficiency  were 
made.  Since  the  model  is  one  dimensional,  the  depth  variations  in  the 
longshore  direction  are  averaged  out.  Sand  bars  and  other  longshore  features 
are    taken    into    account    but    offshore    trenches    and    channels    are    smoothed   out 
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by    this    process.       The    new    height    module,     with    this    simplification,    runs    six 
times  faster   than   the   original   module. 

The  new  wave  height  subroutine  has  many  advantages.  It  uses  a  simple 
but  robust  probabilistic  approach  to  the  wave  height  problem.  It  allows 
variable  bathymetry  in  the  offshore  direction  and  obliquely  impinging  wave 
fronts.  It  is  fast,  accurate,  and  the  governing  equations  can  easily  be  fine 
tuned  by'  varying  the  parameters. 
D.      LONGSHORE    CURRENT    CALCULATION 

The  original  longshore  current  model  (Chen  and  Wang,  1983)  used  a  two 
dimensional  momentum  flux  balance  which  was  numerically  intensive.  The 
model  converged  very  slowly  (over  23  hours  of  CPU  time  in  some  cases),  and 
the  calculations  exhibited  unstable  behavior  (Devendorf,  1985).  To  provide 
the  necessary  longshore  current  information  to  an  operational  planner,  a 
simple  one  dimensional  model  was  substituted  for  the  original  two  dimensional 
subroutine .  The  model  was  designed  to  work  with  the  one  dimensional  wave 
height  model  discussed  in  the  previous  section .  The  random  wave  probability 
model  is  used  as  input  for  the  longshore  current  calculations.  A  narrow 
banded  wave  distribution  is  assumed  so  that  the  wave  frequencies  can  be 
approximated  by  a  mean  frequency  f  and  direction  6  .  The  longshore 
momentum  equation  is  a  balance  between  the  mean  wave  momentum  flux 
divergence  and  the  longshore  component  of  the  bed  shear  stress  (Thornton 
and    Guza,    1985): 


sx  y 
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where : 


>     =   time  and  depth  averaged  covariance 

between  unsteady  velocity  components 

b 

=   bed  shear  stress  in  longshore  direction 


S  is     made     up     of     two     terms     which     are     assumed     to     be     statistically 

yx 

independent.       The   first   term    is    the    radiation   stress    term    (S       )    which   is   the 
wave     induced     momentum     flux.  The     second     term      (S      ')      is     the     depth 

y* 

integrated    turbulent     Reynold's    stress,     often     parameterized     in    terms    of    an 

eddy  viscosity   term. 

Monochromatic    current    models    require    the    eddy    viscosity    term    to   smooth 

out    the    discontinuity     in    the    current    field    at    the    breaker    line.        Since    no 

waves    break    outside    the    surf   zone    in    a    monochromatic    model,     there    is    no 

current    generated.       However,    there    is    a    jump    discontinuity    at    the    breaker 

line    which    requires    smoothing    by    the    eddy    viscosity    term.       In    the    random 

wave    model,    there    is    no    clearly    defined    breaker    line.        Some    waves    break 

offshore    and,    as    the    beach    is    approached,     more    waves    break    until,    in    the 

inner    surf   zone,    almost    all    waves    are    breaking    (Thornton    and    Guza,    1985). 

The    smooth    transition    of   energy    dissipation    reduces    the    requirement    for    the 

eddy  viscosity   term,    which  simplifies  the  velocity   calculation. 

The     radiation    shear    stress     S  is    described    by     (Longuet-Higgins     and 

y  x 

Stewart,    1964): 


S„„      =      EC      cos  6 


yX  **-  ~g    v-uo  v  C  (30) 
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Snell's  law   requires   that: 


—                                               sin    6 
sin   6                     ^                             o  , 

— g —      =      constant      = (31) 

o 


so   that   when    (30)    is  differentiated   the   result   is: 

j       -  sin    9       ,  sin    a" 

dX,  byx  C  dX(ECgcos9)        =       c <Er/         (32) 


This    implies    that    the    shear    stress    is    balanced    by    the    bore    dissipation    for 
breaking   waves   on   a  straight   and  parallel  beach. 
The  bed  shear  stress  is  given   as: 


Ty      =       P  Cf  |5|    V  (33) 


where: 

|U|   =   instantaneous  total  bottom  velocity 
V    =   longshore  current  velocity 

The   bed    stress   is    linearized    by    assuming    small    angles    of    wave    incidence    and 
weak  longshore   currents    (Thornton,    1970;    Longuet-Higgins,    1970): 


y        =         P    Cf   |U|      V  (34) 


Using  shallow   water   wave   theory,    the   bottom   velocity   is   computed  by: 
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Making  the   appropriate   substitutions   into    the   longshore   momentum   relation   and 
solving  for   V  : 

,  c3  -  1/2      -r   H6 

3    B  f  g     sin  6   rms 

=   4    Cf/      Co   h9/2  (36) 


The  numerical  scheme  is  compared  to  the  analytical  solution,  which  is 
obtained   *by    simply    substituting    the    expression    for    H  (Equation    27)    into 

the  velocity  calculation  (Equation  36).  The  numerical  and  the  analytical 
solutions  for  the  shallow  and  steep  wave  cases  studied  above  are  compared  in 
Figures  3.25  and  3.26.  The  numerical  solution  of  V  slightly  underpredicts 
the    current    velocity    compared    with    the    analytic    solution.        This    is    due    in 

part    to    the   introduction    of  a   small   error   in    H  which    was   discussed   in    the 

r  rms 

previous   section.      The  value  of  H  is  raised  to  the  sixth  power  in  Equation 

(36)      which     means     that     V     may     be     sensitive     to     small     errors     in     H 

J  rms 

However,  the  magnitude  of  the  error  for  V  is  on  the  order  of  0.05  and  0.07 
m/sec   which   may   be  neglected. 

The     initial     runs     of     the     updated     wave     height     and     longshore     current 
subroutines    are    presented    and    the    reworked    SURFCON    module   is  discussed   in 
the  next  section . 
E.      IMPROVEMENTS   TO    THE    SURFCON   MODULE 

During  the  reprogramming  effort,  a  major  consideration  was  maintaining 
the  integrity  of  the  original  model,  while  incorporating  modular 
improvements.  The  new  height,  direction  and  current  subroutines  were  coded 
under     new      names,      leaving      the      old      subroutines      intact.  To      maintain 

consistency,  the  control  module,  SURFCON,  was  renamed  OLDSURFCON  and 
a     new     module,      NEWSURFCON,      was     written     to     make     use     of     the     new 
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subroutines.  The  advantage  of  this  approach  is  that  the  model  can  be 
changed  from  the  new  to  the  old  version  simply  by  changing  the  main 
program    "CALL    NEWSURFCON"    to    "CALL    OLDSURFCON". 

The  NEWSURFCON  module  is  basically  the  same  as  the  old  module  with 
some  exceptions.  The  wavenumber  k,  celerity  C  and  group  speed  Cg  are 
initially  calculated.  Then  the  new  wave  direction  and  height  subroutines  are 
called,    which   return   direction,    height   and   longshore   current   matrices. 

The  first  major  change  incorporated  in  this  new  module  is  the  way  the 
breaker  line  is  defined.  In  the  old  module,  a  quadratic  equation  was  solved 
for  the  location  of  the  first  breaker  line.  Since  the  new  breaker  heights  are 
calculated  using  a  weighted  probability  density  function,  Ph(H)  =  W(H)  p(H), 
there  is,  by  definition,  no  specific  spot  where  all  waves  are  said  to  be 
breaking.  The  weighting  function,  W(H),  as  defined  in  Equation  (37),  gives 
a  measure   of  the   percent   of  breaking   waves    (Thornton   and    Guza ,    1983): 


W(H)       =        [%£]  =       /       Pb(H)d(H)  (37) 


To  specify  a  "breaker  line" ,  the  model  first  differentiates  between  steep 
waves,  which  show  no  wave  growth  due  to  shoaling  towards  the  shore,  and 
shallow  waves,  which  increase  in  wave  height  until  they  break.  The  critical 
wave  steepness,  determined  from  Equation  (27),  is  H  /L  =0.02.  A  wave 
with  a  steepness  greater  than  this  criteria  is  considered  steep  while  a  wave 
with  steepness  less  than  0.02  is  considered  shallow.  For  steep  waves,  a 
decision  must  be  made  as  to  what  percentage  of  waves  having  broken 
constitutes  the  surf  line.  After  some  rather  ad  hoc  experimentation,  a  value 
of    0.33    was    chosen    as    a    reasonable,    which    in    effect    defines    a    "significant 

77 


breaker  height".  That  is,  once  33  percent  of  the  waves  are  breaking,  the 
model  chooses  that  point  along  the  grid  as  the  first  breaker  line.  This  test 
is  actually  conducted  in  the  wave  height  subroutine.  As  the  model  moves 
down  the  column,  the  parameter  W(H)  is  tested  at  each  grid  point  to  see  if 
it  greater  than  or  equal  to  0.33.  Once  this  value  is  reached,  the  breaker 
height,  breaker  depth,  and  distance  from  shore  are  returned  to  the 
NEWSURPCON  module.  For  shallow  waves,  the  model  selects  the  maximum 
wave  height  as  the  breaker  line.  These  two  methods  give  good  estimates  of 
the   observed   breaker   line,    as  discussed   in    Chapter   5. 

Once  the  initial  breaker  height  is  defined  in  this  manner,  the  number  of 
additional  lines  are  computed  by  dividing  the  initial  breaker  distance  by  the 
wave  length  corresponding  to  the  average  surf  zone  depth.  If  there  are  one 
or  more  wave  lengths  between  the  initial  breaker  line  and  the  beach,  the 
number      of     breaker      lines      are      incremented      accordingly.  No      additional 

calculations  are  performed  on  these  interior  breaker  lines ,  if  present .  The 
breaker  type  is  calculated  based  on  the  value  of  the  surf  parameter,  Eta, 
which  is  a  function  of  the  beach  slope  and  the  deep  water  wave  steepness 
(Battjes,    1974): 


c,       = 


tan    6  tan 


(H     /I      )1/2  Hn         1/2  (38) 

(Ho/LoJ  (2tt( — ^-))1/2 

gT 


If  Eta  is  less  than  or  equal  to  0.4,  the  breaker  type  is  defined  as  spilling. 
Eta  values  between  0.4  and  2.0  are  defined  as  plunging  and  values  greater 
than  2.0   are   called  collapsing  or  surging   waves. 
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The    next    calculation    performed    by    NEWSURFCON    is    the    "Effective    surf" 

computation.        Effective    surf    is    a    measure    of   surf    intensity    and    is    used    by 

operational    planners    to    provide    a    criterion    for    the    feasibility    of    conducting 

amphibious    landing    operations    using    specific    kinds    of    equipment.        Effective 

surf,     expressed    as    a    wave    height    in    feet,     is    a    modification    of    significant 

breaker  height,    taking  into  account   the   total  beach   and   surf  conditions.    It   is 

important  Xo   note    that    effective   surf  is   simply    a   planning   parameter,    and  not 

a    direct    correspondence    to    the    significant    breaker    height.       Most    commonly, 

H  f    is    larger    than    H   ,    but    in    some    cases,    especially    when    the    H      is    small 

and   the  breaking   waves  are  all  spilling,    H   „  can   be  smaller   than   H   . 

eff  s 

The  calculation  of  effective  surf  uses  a  simple  "look-up  table"  method 
which  is  outlined  in  the  Joint  Surf  Manual  ( COMNA  VSURFPA  C/ 
COMNAVSURFLANT  Inst.  3840.1,  1976).  The  calculations  of  effective  surf 
are  retained  from  the  original  model  with  only  minor  modifications. 
Originally,  the  effective  surf  was  retained  as  a  wave  height  matrix. 
However,  H  „  is  simply  a  single  number  for  a  specific  beach,  based  on 
significant  wave  height  and  direction ,  period ,  wind  speed  and  direction ,  and 
percentage  of  spilling  waves.  The  new  module  calculates  a  single  number  for 
H  ff  and  outputs  all  the  H  „  parameters  as  specified  by  the  Joint  Surf 
Manual.  Other  information  of  use  to  the  planner  is  printed  at  the  same 
time.  The     new     calculation     also     allows     for     a     wind     input      (speed     and 

direction)    which   was  not   available  previously. 

The  final  operation  in   NEWSURFCON   is   the   limitation   of  the   wave   heights 
to    the    maximum    stable    wave    height    for    a    given    depth    and    period.        This 
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maximum   wave   height   is  given   as    (Wang   and    Chen,    1983): 


H      Y      =      0.067  tanh  (kh)  (^)  (39) 


Each    generated    wave    height    is    compared    against    this    maximum    value    and    is 

reset    to    H  if    it    exceeds    the    limiting    wave    height.       The    wave    height, 

max 

direction,    current    and   effective    surf   matrices    are    then    sent    to    OUTSURF    for 
formatted   output   to   either   the   screen   or   to   a   designated  printer. 
F.      DEVELOPMENT    OF    AN    OPERATIONAL    USERS    GUIDE 

A  simple  users  guide  was  written  for  the  new  version  of  the  SSSP  (Gill, 
1985).  Several  guidelines  were  adhered  to  during  its  development.  The 
manual  assumes  some  familiarity  with  the  HP-9845/275  and  a  working 
knowledge  of  oceanographic  terms.  The  manual  is  tutorial  in  nature.  It  is 
divided  into  two  parts,  for  the  two  main  modules  OPEN  and  SURFCON.  The 
manual  explains  each  input  request  and  appropriate  examples  are  included 
where  necessary.  An  effort  was  made  to  keep  the  manual  brief  and  user 
friendly.  The  model,  being  extensively  menu  driven,  is  straight-forward 
enough  to  require  little  documentation,  after  it  has  been  used  once  or  twice. 
The  BASIC  code  is  extensively  documented  throughout.  Where  appropriate, 
the  old  program  lines  are  left  in  place  (commented  to  prevent  execution), 
with  the  replacement  lines  documented  to  point  out  the  reason  for  the 
specific  change. 
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IV.      DISCUSSION 

A.      SENSITIVITY   TEST    OF    BREAKER    PARAMETERS 

As  discussed  in  Chapter  III,  the  computation  of  H  ,  using  the  Rayleigh 
probability  density  function  model,  depends  on  two  breaker  parameters, 
Gamma  and  B.  Gamma  is  an  adjustable  coefficient  of  0(1)  which  appears  to 
be  strongly  related  to  the  beach  slope  (Sallenger  and  Holman,  1985).  B  is 
the  breaker  coefficient  which  represents  the  percentage  of  foam  on  the  face 
of  the  breaking  wave  (i.e.,  the  intensity  of  the  breaking  wave).  This 
parameter  is  expected  to  be  less  than  or  equal  to  1.0,  since  a  B  value  of 
1.0  corresponds  to  a  fully  developed  bore  (a  wave  that  has  fully  broken  from 
crest  to  trough).  Values  of  B  greater  than  1.0  imply  that  the  bore 
dissipation  function  underestimates  the  wave  dissipation  due  to  breaking 
(Thornton   and   Guza ,    1983). 

To  verify  the  validity  of  the  model  assumptions,  24  wave  data  sets  were 
studied  to  gain  an  understanding  of  the  interaction  of  the  B  and  Gamma 
values  as  used  in  the  model.  These  data  sets  are  summarized  in  Table  4.1. 
The   data    were   digitized   and   appropriate    Gamma   values   chosen,    based    on    the 

slope   of  H        /h      inside    the   surf   zone.      The    Gamma   values   ranged   from    0.39 

r  rms 

to  0.60  and,  for  each  data  set,  the  specific  Gamma  value  was  held  constant 
while  the  B  values  were  varied.  A  curve  fitting  routine  was  used  to 
optimize  the  B  parameter  for  a  given  Gamma.  As  B  was  varied,  the  new 
model  curve  was  compared  with  the  real  data  points,  and  a  sum  of  the 
square  of  the  errors  was  generated.  The  B  values  were  iterated  until  the 
error  sum  was  minimized.  Plots  of  the  model  output  versus  the  data  points 
were  generated  for  each  case. 
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A    total   of  24    data   sets   were    tested   against   the   model    (See  Table   4.1). 

In    every    case,    the    model   accurately    predicts    the   wave   heights  described   by 

the    data    points.       The    first    ten    data    sets    are    from    a    laboratory    study    of 

energy    saturation    of    irregular     waves     during    shoaling     (Vincent,     in    press). 

These    wave    tank    tests    were    conducted    with    a    plane    bottom    slope    of   1:30, 

with    various   frequency    and    incident    wave    heights.       The    experiments    can    be 

divided    up    into    two    groups    based    on    wave    steepness.       The    low    steepness 

cases     (H    /L       <     0.02)     are     represented    by     Case     1638     (Fig.     4.1).        The 
o     o 

Gamma    value    of    0.60    was    derived    from    the    slope    of   H  /h    as    described 

c  rms 

above.  Using  the  curve  fitting  routine,  an  optimum  B  value  of  0.82  was 
chosen.  Note  the  good  fit  of  the  model  curve  (solid)  compared  with  the 
data  points  (circles).  The  high  steepness  cases  are  represented  by  Case 
1148  (Fig.  4.2).  Gamma  was  set  equal  to  0.60  and  a  B  value  of  1.2  was 
obtained.  As  before,  there  is  good  agreement  between  the  model  curve  and 
the  data  points. 

To  determine  the  optimized  model's  behavior  when  used  with  more 
difficult  bathymetry,  data  from  several  wave  tank  experiments  with  offshore 
bars  were  tested  (Battjes  and  Stive,  in  press).  Case  5  (Fig.  4.3)  has  a 
simple  offshore  bar  with  a  mean  bottom  slope  of  1:25.  Gamma  was  estimated 
to  be  .45  and  a  B  value  of  1.2  was  obtained.  Although  the  presence  of  the 
offshore  bar  (Fig.  4.4)  radically  changes  the  wave  height  data  profile,  the 
optimized  model  curve  again  faithfully  predicts  the  wave  heights.  Wave  tank 
test  case  15  (Fig.  4.5)  has  a  very  complex  bathymetry  (Fig.  4.6).  The 
model  divides  the  on /offshore  distance  into  10  grid  points  with  the 
intermediate    depths    calculated    by    straight    linear    interpolation.       Therefore, 
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complex  bathymetry,  such  as  that  in  Case  15,  was  smoothed  out  to  a 
certain  degree.  Even  with  that  limitation  in  mind,  the  optimized  model 
curve  appears  to  accurately  estimate  the  recorded  wave  heights.  The  model 
results   will  be   tested  against  field  data  in  the  next  chapter. 

It  is  evident  that  the  random  wave  model,  with  the  appropriate  choice 
of  the  surf  parameters,  Gamma  and  B,  will  accurately  forecast  wave  heights 
throughout-  the  domain.  A  sensitivity  test  for  the  choice  of  parameters  was 
conducted  on  three  representative  cases  1638,  1148,  and  5.  Case  5  (Battjes 
and  Stive,  in  press)  was  chosen  for  discussion  here  because  of  its  more 
interesting     bottom     profile.  It     was     felt     that     this     barred     profile     was 

representative  of  many  real  beaches  and  would  provide  a  more  rigorous  test 
than  a  simple  plane  bottom.  Using  a  Gamma  value  of  0.45,  the  optimum  B 
value  of  1.3  was  changed  to  1.7,  1.5,  1.1  and  0.9  (Fig.  4.7).  B  values 
less  than  optimal  tend  to  over-estimate  the  wave  heights  and  vice-versa.  It 
is  of  interest  to  note  that  changes  in  B  of  22-33  percent  only  results  in  an 
increase  in  the  model  error  of  about  10-12%.  This  was  confirmed  in  a  study 
of  Torrey  Pines,  California  wave  data  by  Thornton  and  Guza  (1983).  They 
found  that  a  variation  of  +/-  25%  about  the  optimum  B  value  resulted  in  an 
increased  model  error  of  less  than   10%   (Fig.    4.8). 

It  was  found  that  the  model  demonstrated  a  similar  sensitivity  to  changes 
in  the  Gamma  value.  A  change  of  about  +/-  20%  to  the  optimum  Gamma 
value  of  0.45  resulted  in  an  increased  model  error  of  8-10%  (Fig.  4.9). 
Sensitivity  tests  were  conducted  on  high  and  low  steepness  plane  beach  cases 
(1638  and  1148)  with  similar  results.  In  the  next  section,  techniques  and 
rationale  for  choosing  B  and  Gamma  values,  based  on  initial  conditions,  are 
discussed.       It    is    important    to    bear    in    mind    the    above    sensitivity    analysis. 
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The    model    appears    to    be    fairly    robust    with    respect    to    small    changes    in    the 
surf    parameters.         Therefore,     slight     errors     in     the     choice     of    these    values 
should  not   degrade   the  final  performance   of  the   model   to   a   great   extent. 
B.       CHOOSING    THE    BREAKER    PARAMETERS 

The    two    breaker   parameters,    B    and    Gamma,    are    multiplicative    constants 

in     the     H  wave     height     calculation      (see     Equation     27).         They     could 

rms  J 

therefore t  be  combined  as  a  single  bulk  modulus  coefficient,  or  kept  as  two 
separate  parameters.  It  was  felt  that  keeping  the  two  terms  separate  would 
better  serve   the   requirements  of  the   model. 

Some  work  has  already  been  done  on  the  problem  of  choosing  the  breaker 
parameters.  A  recent  study  by  Sallenger  and  Holman  (1985)  has  determined 
a  functional  relationship  between  Gamma  and  the  average  beach  slope  based 
solely  on  field  observations  inside  the  surf  zone .  This  relationship  is  given 
by: 


Y      =      3.2  tan  6    +    0.3 

(40) 


Using  the  24  wave  data  sets,  the  chosen  Gamma  values  are  plotted  against 
beach  slope  (Fig.  4.10).  The  solid  line  is  a  straight  linear  regression  of  the 
data  points  while  the  dashed  line  is  the  Holman  /Sallenger  function.  There 
seems  to  be  a  general  agreement  in  the  slope  of  the  line  although  the 
intercept  is  slightly  different.  Gamma  is  also  plotted  against  wave  steepness 
(Fig.  4.11)  but  the  plot  is  less  linearly  correlated.  As  an  initial  trial,  the 
Holman /Sallenger  functional  relationship  will  be  adopted  as  the  model's 
convention  for  choosing   the    Gamma   value.      The   average  slope  is  defined  here 
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as  the  depth  at  the  outermost  grid  point  divided  by  the  offshore  grid 
distance . 

Choosing  a  proper  predictor  for  the  B  value  proved  to  be  more  difficult. 
Since  B  represents  a  measure  of  breaker  intensity,  it  is  felt  that  it  would  be 
related  to  incident  wave  steepness  which  has  previously  been  shown  to  be 
related  to  breaker  type  (Battjes,  1972).  The  optimum  B  values  for  the  24 
wave  sets  are  plotted  against  steepness  in  Figure  4.12  and  against  beach 
slope  in  Figure  4.13.  The  wave  steepness  in  this  case  was  defined  as  the 
incident  wave  height  divided  by  the  deep  water  wave  length  (H  /L  ). 
Although  there  is  a  considerable  scatter,  B  appears  correlated  with  wave 
steepness  and   less   with   beach  slope. 

Several    other    parameters    were    considered    as    predictors    for    B .       The   B 

and    Gamma    values    were    compared    with    the    surf   parameter,    Eta    =    tan    8    / 

1/2  3  4 

(H    /L    )         .         Additionally,     the     Bulk     Modulus     (B    /Gamma    )      was     plotted 

against   beach  slope   and   wave   steepness.      Since,   from    Figures  4.11    and   4.12, 

B    appears    to    be    related    to    both    beach    slope    and    wave    steepness,    it    was 

hoped   that   Eta,    which   is  a  measure   of  the    expected   breaker   intensity,    would 

offer    a    predictor    for    B.        The    plot    of    B    versus    Eta    (Fig.    4.14)    shows    an 

interesting    pattern.       B    appears    to    increase    almost    linearly    as    Eta    increases 

to  about   0.6   and  then   decreases  in   an   exponential  fashion   to   an   asymptotic   B 

value    of    0.8.        This    behavior    appears    reasonable.        As    the    surf    parameter 

increases,    indicating   a   transition  from    spilling    to   plunging   waves,    the   B   value 

(amount   of  wave  face   that   is  breaking),    should  increase. 

The    maximum    value    that    B    should    reach,     according    to    theory,    is    1.0. 

Studies    have    shown,    however,    that    in    many    cases,    the    optimum    B    value    is 

on    the    order    of    1.4     (Thornton    and    Guza,     1983).        The    maximum    B    value 
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occurs  when  Eta  is  about  0.6,  indicating  that  the  breakers  are  mixed  spilling 
and  plunging.  As  Eta  continues  to  increase,  the  B  values  fall  more  into  line 
with  theory  and  approach  the  expected  value  of  1.0,  as  the  waves  become 
fully   plunging. 

To  simplify  and  expedite  the  model  development,  the  following  piecewise 
continuous  function  is  used   to  determine   the   B   term    (Fig.    4.14): 

B   =     2.167   Eta   +   0.4  Eta   <=  0.6  (41) 

B   =  -   1.70   Eta   +   2.7  0.6   <    Eta   <=  1.0 

B   =     1.0  Eta   >   1.0 

This  choice  is  somewhat  ad  hoc  and  is  based  primarily  on  the  data  sets 
discussed  above.  Future  refinements  in  the  techniques  for  predicting  the 
breaker  parameters  may  substantially  improve  the  use  of  this  model 
operationally.       Using   equations    (40)    and    (41)    as   predictors   for    Gamma   and 

B,  the  SURFCON  module  is  essentially  complete.  The  next  chapter  will 
discuss  testing  the  model  against  two  well-documented  natural  beach  data 
sets,  to  get  a  feeling  for  how  the  improved  model  behaves  against  actual 
data. 

C.  MOTIVATION   FOR   THE    RANDOM    WAVE    MODEL    SUBSTITUTION 

As  discussed  in  Chapters  3  and  4,  the  primary  motivation  for  substituting 
the  one  dimensional  random  wave  model  for  the  two  dimensional  model  was  to 
increase     computational     efficiency     and     accuracy.  The      improved     model 

calculates  longshore  current  information  quickly  and  accurately.  Significant 
wave  height  predictions  compare  well  with  actual  data  sets,  as  will  be 
discussed  in   Chapter  5.      There   is    ,    however,    a  potential  for   loss  of  detailed 
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information  by  using  the  one  dimensional  model.  Nearshore  circulation  and 
longshore  variation  in  wave  heights  are  smoothed  out  by  the  one  dimensional 
approach . 

The  surf  model  (both  old  and  new  versions)  uses  a  single  wave  height 
input  to  define  the  offshore  boundary  condition.  This  implies  that  variations 
in  the  interior  wave  heights  and  directions  are  due  only  to  variations  in 
bathymetry.  Because  the  model  uses  a  very  small  grid  (only  13  longshore 
grid  points),  very  accurate  bathymetry  must  be  entered.  It  is  questionable 
whether  the  operational  planner  would  have  access  to  bathy  metric  data  of 
the  required  accuracy  to  define   a  two  dimensional  bathymetry. 

Many  beaches  which  would  potentially  be  used  for  amphibious  operations 
are  shallow,  almost  planar  in  nature.  Other  common  landing  beaches  are 
barred,  but  rather  uniform  in  the  longshore  direction.  Therefore,  for  many 
operationally  relevant  beaches ,  the  assumption  of  straight  and  parallel 
contours,  and  the  use  of  a  one  dimensional  model,  appears  justified.  The 
"mean  beach"  surf  information  provided  by  the  model  will  adequately 
represent  the  conditions  over  a  standard  500  yard  landing  beach.  If  a  wider 
beach  needs  to  be  characterized,  the  improved  model  would  be  run  several 
times ,    entering   the   required  information  for  each   beach  segment . 
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V.      TEST  AND   EVALUATION   OF  THE  IMPROVED   MODEL 

A.       WAVE    HEIGHT    AND    BREAKER    LOCATION 

The  improved  version  of  the  SSSP  was  tested  against  one  wave  tank 
experiment  and  two  sets  of  beach  surf  data  to  see  if  the  model  would  choose 
reasonable ;  Gamma  and  B  parameter  values,  based  on  the  predictors  discussed 
in  Chapter  4.  The  predicted  wave  heights  and  surf  characteristics  are 
compared   with   the  observed  values. 

Case  5  (Battjes  and  Stive,  1985),  discussed  in  the  preceding  chapter, 
was  chosen  as  the  wave  tank  experiment  with  which  to  test  the  finished 
model.  This  case  was  chosen  because  of  its  barred  beach  bathymetry. 
The  model  selected  values  of  0.41  and  0.963  for  Gamma  and  B.  The 
optimum  values,  from  Table  4.1,  are  0.45  and  1.30.  The  model  wave  height 
fit  to  the  observed  data  points  is  still  quite  reasonable  (Fig.  5.1).  The 
solid  line   represents  the  model  output  and  the  circles  are   the  observed  data. 

Two  Southern  California  beach  data  sets  were  used  for  further  tests. 
Data  collected  at  Santa  Barbara's  Leadbetter  Beach  during  February,  1980 
and  Torrey  Pines  during  November,  1978  (Thornton  and  Guza,  1983),  were 
used  as  the  primary   tests  of  the  performance  of  SSSP. 

The  first  Santa  Barbara  test  used  data  collected  on  3  February,  1980. 
The  average  beach  slope  was  0.047  and  the  incident  waves  were  0.55  meters 
(rms)  with  a  period  of  14.3  seconds.  The  observed  breaker  height  was 
measured  at  0.70  meters,  approximately  40  meters  from  the  beach  (dashed 
line).      The   wave   angle   was   7.8   degrees  from   normal.      The  SSSP  predicted  a 
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rms  breaker  height  of  0.68  meters  (solid  line),  40.0  meters  from  the  beach 
(Fig.  5.2).  The  model  under  predicted  the  observed  heights  less  than  4 
percent . 

Similar  results  are  obtained  for  Leadbetter  Beach  data  obtained  on  4  and 
5  February,  1980.  On  4  February,  the  incident  wave  height  was  0.56 
meters  with  a  period  of  14.3  seconds.  The  wave  angle  was  9  degrees  from 
normal.  The  observed  breaker  height  was  0.67  meters,  50  meters  from  the 
beach.  The  model  predicted  0.65  meters,  40  meters  from  the  shore  (Fig. 
5.3). 

On  5  February,  the  incident  wave  was  0.45  meters  with  a  period  of 
12.8  seconds.  The  wave  angle  was  8.4  degrees  from  normal.  Observed 
breakers  were  0.60  meters,  about  36  meters  from  the  beach.  SSSP  predicted 
0.55    meter  breakers,    33   meters  from   the  shore    (Fig.    5.4). 

Torrey  Pines  beach  data  acquired  on  4  and  10  November,  1978  is 
compared  with  model  results  in  Figures  5.5  and  5.6.  The  beach  has  a  shallow 
slope  of  0.022  which  accounted  for  a  wide  surf  zone  of  about  160  meters. 
On  4  November,  the  incident  waves  were  about  0.35  meters  with  a  period  of 
14.3  seconds  and  normal  incidence.  The  observed  breaker  height  was  0.47 
meters  about  110  meters  from  the  shore.  SSSP  predicted  0.48  meter 
breakers  at  a  distance  of  106  meters  from  the  beach.  The  model  correctly 
forecast  the  wave  height  and  surf  zone   width   (Fig.    5.5). 

A  similar  result  is  seen  for  the  November  10  Torrey  Pines  Case.  The 
incident  wave  is  0.68  meters  with  a  period  of  15.9  seconds  and  normal 
incidence.  Observed  breakers  were  0.72  meters,  approximately  130  meters 
from  the  shoreline.  The  model  predicted  breaker  heights  of  0.75  meters,  140 
meters   from    the   beach    (Fig.    5.6).      As  for   the   November   4    case,    the   model 
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worked  quite   well  for   wave   height   prediction   and   the   location   of  the  breaker 
line. 

Based  on  the  tests  of  the  SSSP  against  various  data  sets,  it  appears  that 
the  new  version  is  an  improvement  over  the  original  model.  The  old  version 
consistently  underpredicted  the  wave  heights  and  was  very  inaccurate  at 
estimating  the  location  of  the  breaker  line  (Devendorf,  1985).  The  new 
version  appears  to  be  quite  good  at  predicting  the  breaker  height  and  gives 
better  estimates  of  the  breaker  line.  Table  5.1  is  a  comparison  between  the 
optimum  values  of  the  two  breaker  parameters  and  the  values  chosen  by  the 
model.  The  technique  for  choosing  Gamma  and  B  is  convenient,  although 
improvements  will  need  to  be  made  in  this  area. 
B.      LONGSHORE    CURRENT    DISTRIBUTION 

Longshore  currents  are  the  result  of  waves  impinging  at  an  oblique  angle 
to  the  shoreline.  All  the  wave  tank  tests,  as  well  as  the  Torrey  Pines 
cases,  had  normally  incident  waves.  The  Santa  Barbara  Leadbetter  Beach 
data  is  the  only  comprehensive  field  data  on  longshore  currents  acquired  to 
date.  There  exists  no  laboratory  data  for  longshore  currents  generated  by 
obliquely  incident,  random  waves.  The  same  three  data  sets  (February 
3,4,5)  were  chosen  as  study  cases  because  they  had  a  variety  of  incident 
wave  heights  and  wave  angles.  For  the  February  3  case  (Fig.  5.7),  the 
maximum  observed  current  (dashed)  was  0.48  m/sec  at  a  distance  of  22 
meters  from  the  shoreline.  The  bed  shear  stress  coefficient,  Cf,  was  set 
equal  to  0.009  for  all  cases.  This  value  was  derived  experimentally  by 
Thornton  and  Guza  (1985).  The  SSSP  prediction  (solid)  for  longshore  current 
was  a  maximum  of  0.5   m/sec,    25   meters  from   the  beach   (solid  line). 
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For  the  February  4  case  (Fig.  5.8),  the  maximum  observed  current  was 
.48  m/sec,  25  meters  from  shore;  the  model  prediction  was  .44  m/sec  at 
20  meters  from  the  beach.  The  Feb  5  case  (Fig.  5.9)  shows  an  observed 
current  maximum  of  .36  m/sec  at  20  meters  from  the  beach.  The  model 
predicted    .36   m/sec  at  20   meters  from   the  shoreline. 

The  model  not  only  predicts  the  maximum  current  velocity  reasonably 
well,  but  the  current  distribution  within  the  surf  zone  shows  good  agreement 
with  the  data.  The  model  seems  to  underpredict  the  current  velocity 
seaward  of  the  maximum  current,  but,  in  general,  the  model  does  quite  a 
good  job  of  longshore  current  prediction. 

The  other  Leadbetter  Beach  cases  showed  similar  results.  In  all  cases, 
the  model  gives  reasonable  maximum  current  locations  and  speeds  as  well  as 
current  distributions  within  the  surf  zones.  These  findings  were  also 
confirmed  by  Thornton  and  Guza  (1985)  in  a- study  in  which  the  random  wave 
model's  wave  height  and  longshore  current  predictions  are  rigorously  tested 
against   the   Leadbetter  Beach  data. 
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VI.      SUMMARY   AND   CONCLUSIONS 

A.      SUMMARY 

In  1981,  the  U.S.  Navy  contracted  to  have  a  Sea,  Swell  and  Surf 
Program  (SSSP)  developed  to  provide  input  to  the  planning  efforts  of  Naval 
Oceanographers  at  sea  and  amphibious  operational  planners.  The  model  was 
developed  in  FORTRAN  and  designed  to  run  on  a  main  frame  computer  (Wang 
and  Chen,  1983).  The  completed  SSSP  was  then  translated  into  BASIC  for 
use  on  the  currently  available  shipboard  micro-computers,  the  HP-9845B/275 
(Devendorf,      1985).  The     SSSP     was     modular     in     nature,      with     separate 

subroutines  for  offshore  wind  wave  generation ,  surf  zone  calculations  and 
nearshore  circulation. 

The  original  SSSP  had  several  problems  which  prevented  its  use  by 
operational  planners.  The  offshore  wind  wave  generation  module  overbuilt 
the  wave  heights  for  high  wind  speeds.  Some  of  the  surf  condition 
calculations  were  inaccurate  and  the  nearshore  circulation  subroutine  had 
stability  problems.  These  problems  were  addressed  by  this  thesis  in  an 
attempt   to  improve  and  upgrade   the   SSSP. 

The  open  water  module  shows  satisfactory  wind  wave  generation  for  wind 
speeds  of  15  knots  or  less.  For  higher  wind  speeds,  the  predicted  wave 
generation   occurs    too   quickly,    both   spatially    and   temporally.  Sensitivity 

analysis  indicated  that  the  Phillips  linear  term  was  responsible  for  the  rapid 
development,  but  no  improvements  to  the  model  physics  were  made.  Several 
ad  hoc  "fixes"  were  recommended.  This  problem  will  require  further 
research . 
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Two  subroutines  are  called  to  calculate  the  wave  transformation  as  the 
wave  proceeds  from  offshore  to  onshore.  Several  numeric  improvements  were 
made  which  eliminated  minor  instabilities  in  the  wave  direction  and  height 
subroutines.  The  original  no  flux  boundary  conditions  were  changed  to  cyclic 
boundary  conditions  in  both  subroutines.  This  change  allows  the  wave 
heights  and  directions  to  more  accurately  reflect  the  conditions  required  by 
the  model  physics  at   the  model  boundaries. 

The  most  serious  problem ,  the  instability  of  the  nearshore  circulation 
module,  is  solved  by  substituting  a  one  dimensional  longshore  current  model 
for  the  original  two  dimensional  scheme.  Although  there  is  a  loss  of  two 
dimensional  current  resolution  within  the  surf  zone  due  to  the  assumptions  of 
straight  and  parallel  contours,  stationary  wave  conditions,  and  shallow  water 
initial  waves,  it  was  felt  that  a  model  that  gives  the  important  longshore 
current  information  quickly  and  accurately  is  an  improvement  over  the 
original  unstable,  two-dimensional  model.  The  new  current  model  is  a 
simple,  longshore  momentum  flux  balance  equation.  It  retains  the  wave 
induced  momentum  flux  and  neglects  the  turbulent  Reynold's  stress.  The  new 
model  is  fast  and  accurate,  requiring  only  a  few  calculations  once  the  wave 
height  field  is  generated. 

To  take  advantage  of  state  of  the  art  breaker  height  modeling,  a  random 
wave  model  was  implemented  for  the  calculations  of  the  wave  heights  in  the 
nearshore  area.  The  random  wave  model  assumes  that  wave  heights,  in  the 
nearshore  region,  follow  a  Rayleigh  distribution  and  that,  even  after  the 
waves  break,  this  distribution  is  retained. As  a  future  project,  the  random 
wave  model  could  be  developed  into  a  two  dimensional  model  by  applying  the 
proper  longshore  numerical  differencing. 
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The  random  wave  model  assumes  that  wave  breaking  can  be  described  by 
a  periodic  bore.  Two  parameters  are  introduced  by  the  model,  which  must 
be  chosen  based  on  the  initial  conditions.  The  Gamma  term  is  part  of  the 
proportionality  coefficient  which  modifies  the  initial  Rayleigh  distribution,  and 
is  found  to  be  strongly  correlated  with  the  beach  slope.  The  breaker 
parameter ,  B ,  describes  the  amount  of  the  wave  face  that  is  breaking  and 
appears  to  be  related  to  the  surf  parameter  Eta.  The  criterion  for  choosing 
the   two  breaker  parameters  was  developed  empirically. 

The  wave  height  and  current  models  were  next  tested  against  wave  tank 
and  natural  beach  data  sets.  Although  the  model  did  not  necessarily  choose 
optimum  values  for  the  breaker  parameters,  the  model  is  fairly  robust,  and 
slight  deviations  in  the  parameter  values  do  not  degrade  the  wave  height 
calculations  significantly.  Good  agreement  was  found  between  the  observed 
wave  heights  and  the  model's  predicted  values.  The  longshore  current,  based 
on  the  model's  predicted  wave  height  field,  were  also  in  close  agreement 
with  observations. 

The  calculation  of  Effective  Surf,  an  amphibious  planning  guideline,  was 
improved.  The  algorithm  was  implemented  in  accordance  with  the  Joint  Surf 
Manual  (1976)  with  respect  to  its  input  requirements  and  its  formatted 
output . 

A   simple  users  manual,    tutorial  in   nature,    was   developed  for  this  version 
of  the   SSSP.      It   is  short,    and  assumes  some   knowledge   of  the   host  computer 
and  basic  oceanographic  terms.      It  is  written   within  the  guidelines  set  forth 
by   the   Navy   for   development   of  computer   models  for   the   Shipboard  Numerical 
Aid  Program    (Brown,    1984). 
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B.       CONCLUSIONS 

The  improved  Sea,  Swell  and  Surf  Program  (SSSP)  can  provide  much 
useful  information  to  the  Naval  Oceanographer  at  sea  or  amphibious  planner. 
It  has  several  problems  which  have  not  been  resolved,  so  its  use  is 
recommended  with  several  caveats.  The  offshore  wind  wave  generation 
module  is  felt  to  be  accurate  for  low  wind  speeds  only.  Useful  information 
on  offshore  wave  growth  can  be  obtained  for  input  wind  speeds  of  15  knots  or 
less.  The  results  obtained  for  higher  wind  speed  inputs  must  be  examined 
carefully  before  use .  The  requirement  for  providing  an  initial  wave  energy 
spectrum  is  probably  impractical  for  shipboard  use.  The  option  of  providing  a 
single  initial  wave  height  and  period  will  be  of  more  value  to  the  operational 
user. 

The  new  random  wave  breaker  height  model  appears  to  work  well.  The 
breaker  heights  and  locations  as  well  as  the  longshore  current  field  are 
accurately  predicted  by  the  model.  Due  to  its  one-dimensional  nature  and 
the  assumption  of  straight  and  parallel  contours,  longshore  wave  height  and 
current  resolution  is  restricted.  However,  the  assumption  that  an  amphibious 
landing  beach  is  uniform  in  the  longshore  direction  is  a  reasonable  one  over  a 
wide  range  of  operational  areas.  The  small  grid  size  of  the  model  supports 
the  use  of  a  one  dimensional  model  to  compute  an  average  set  of  charac- 
teristics for  a  "mean"  beach.  The  speed  and  stability  of  the  new  model 
make  it  a  viable  tool  for  operational  planners.  Recommendations  for  further 
improvement  of  the   SSSP   include: 

1.  Pursue  the  problem  of  the  overbuilding  of  the  offshore   wind   waves. 

2.  Development  of  better  predictors  for  the  breaker  parameters  Gamma 
and  B  by  including  a  wider  data  base.  As  was  seen  in  the 
optimized    tests,     correct    choice    of    these    parameters    will    provide   a 
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very      good      fit      to      observations,      even      with      very      complicated 
bathymetry. 

3.  Explore  the  feasibility  of  using  the  random  wave  model  two 
dimensionally.  This  would  offer  the  advantages  of  this  more  realistic 
approach  to  wave  height  modeling,  while  retaining  the  details  that 
two  dimensional  spatial  resolution  provides. 

The  improved   SSSP   offers   the   Naval    Oceanographer   at   sea  and  amphibious 

operational    planners    a    new    and    useful    tool    that,    when    used    in    conjunction 

with  surf  observations ,    weather  information ,   and  other  available  data ,    can   be 

a  great  aid  in  the  safe  planning  and  execution  of  beach  landing  operations. 
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